Contraintes

Une contrainte impose une relation aux inconnues du problème (une valeur imposée, une liaison) sans être une physique au sens d’une loi de comportement : elle n’a ni matériau, ni intégrande volumique. pyrucast les traite comme des SubModel ordinaires, contribuant leurs blocs à la même matrice globale que les physiques — pas de système point-selle séparé à orchestrer côté utilisateur.

Ce qu’elles ont en commun

Multiplicateurs de Lagrange

Les contraintes sont imposées par multiplicateurs de Lagrange. Une relation C·u = u_d introduit une inconnue supplémentaire λ (le multiplicateur) et deux blocs rectangulaires dans la matrice :

\[ \begin{bmatrix} K & C^\top \\ C & 0 \end{bmatrix} \begin{bmatrix} u \\ \lambda \end{bmatrix} = \begin{bmatrix} f \\ u_d \end{bmatrix} \]

• le bloc C porte la relation de contrainte (lignes = duale propre, colonnes = primale contrainte de la cible) ;
• le bloc Cᵀ réinjecte l’effort du multiplicateur dans l’équation de la physique cible ;
• à la solution, le multiplicateur λ est l’effort de la pièce sur la contrainte : la première ligne donne Cᵀλ = f − K·u.

La réaction est la force que la contrainte exerce sur la pièce : le champ −Cᵀλ, aux nœuds contraints. Elle se lit par −internal_forces(modele.filter("constraint"), …), que solve ait gardé les multiplicateurs ou éliminé les relations : internal_forces d’une contrainte rend Cᵀλ. Le pas à pas la range dans ses résultats, sous "reactions".

Chacun des blocs C et Cᵀ est, pris isolément, non symétrique ; seule leur union C ∪ Cᵀ l’est — c’est une propriété globale du système point-selle, pas de chaque bloc.

C’est exactement ce que la déclaration de symétrie rend exprimable : chacun des deux déclare porter une moitié (Symmetry::Half), avec une identité de paire commune que seul constraint_block_pair mint — celui qui écrit le même coefficient des deux côtés. L’agrégat conclut à la symétrie quand il tient les deux moitiés, et la retire si un découpage en sépare une.

Les nœuds-multiplicateurs sont des nœuds comme les autres, créés par la contrainte : un par relation, dans le même Coords, à la position d’un nœud de la relation (le nœud imposé pour Dirichlet, celui du premier terme pour une MPC, le nœud immergé pour le baignage, le nœud esclave pour le contact). multiplier_mesh() les donne.

Imposition par élimination (condensation)

Voie alternative aux multiplicateurs, sélectionnée au moment de la résolution par solve(…, relations="eliminate"). Plutôt que d’agrandir le système, on élimine chaque relation : un terme esclave s est exprimé par les maîtres,

\[ u_s = \frac{1}{a_s}\Big(g - \sum_{k \neq s} a_k\, u_k\Big), \]

d’où une transformation globale u = T·û + u₀ (û = DOFs retenus). Le système se réduit à K̂ û = f̂ avec K̂ = Tᵀ K T, f̂ = Tᵀ(f − K·u₀), résolu par la même factorisation, LU ou Cholesky, mais sur une matrice plus petite et définie, sans inconnue multiplicateur. La solution est prolongée u = T·û + u₀, et le multiplicateur de chaque relation est récupéré en post-traitement, aₛ·λ = f − K·u à la ligne de son esclave, puis rendu à son nœud multiplicateur : la solution est celle de la voie Lagrange, mêmes inconnues sur les mêmes supports.

Les deux voies lisent la même description méthode-neutre des contraintes (Constraint::relations()) et acceptent le même chargement (le second membre g au nœud-multiplicateur, ci-dessous) : passer de l’une à l’autre ne change que l’appel de résolution.

Périmètre : non chaîné, esclaves disjoints — chaque relation élimine un esclave distinct, jamais réutilisé comme maître ni esclave dans une autre relation (couvre la périodicité). Un système chaîné est refusé avec une erreur explicite ; il reste résoluble par la voie Lagrange (solve).

Relations unilatérales (inégalités)

Toute relation d’égalité peut devenir unilatérale : Σₖ aₖ·uₖ ≥ g (ou ≤ g) au lieu de = g. C’est le paramètre optionnel sense des contraintes ("=" par défaut, ">=", "<=") — une butée u ≥ a est un Dirichlet unilatéral, une liaison à jeu est une MPC unilatérale :

barre = pyrucast.model.truss(fes)
butee = pyrucast.model.dirichlet(imposed, barre, "u_x", sense=">=")

Une relation unilatérale obéit aux conditions de complémentarité (KKT) : ou bien elle est active (l’égalité tient, et sa réaction pousse), ou bien elle est inactive (le jeu C·u − g est du côté admissible et λ = 0). Avec la convention des blocs assemblés (K·u + Cᵀ·λ = f), le signe admissible du multiplicateur est λ ≤ 0 pour une relation ≥ active, λ ≥ 0 pour une ≤ active.

Comme l’ensemble actif n’est pas connu d’avance, la résolution est itérative — la méthode du statut (active-set), que solve lance dès que la matrice porte une inégalité :

1. statut initial : toutes les inégalités actives (ou le statut convergé précédent quand le cache est chaud — warm start) ;
2. résolution du système point-selle avec, pour chaque relation inactive, sa ligne de contrainte remplacée par λ = 0 (la matrice garde sa taille, seules les valeurs changent) ;
3. mise à jour du statut : une relation active dont le λ tire (signe inadmissible) est relâchée, une inactive dont le jeu pénètre est activée ;
4. statut stable ⇒ convergé ; sinon on refactorise et on répète (boucle finie de la méthode du statut classique, bornée par max_iter).

# La matrice sait quelles relations sont des inégalités : `solve` lance l'ensemble actif.
solution = pyrucast.solver.solve(k, rhs)  # active_set, max_iter, tol

Les relations inactives sortent avec λ = 0 exact ; la solution a la même forme que la voie Lagrange (primal + multiplicateurs). Les coefficients assemblés sont inchangés : le sense est posé sur le bloc C de la contrainte, et c’est là que solve le lit — aucun modèle n’est passé à la résolution, aucun ne peut contredire la matrice. Une matrice sans inégalité est résolue par la factorisation ordinaire.

Le second membre g (la borne) suit le mécanisme habituel : au slot imposed du nœud-multiplicateur, donc constraint_rhs fonctionne tel quel.

Second membre : le helper constraint_rhs

Le second membre u_d / g n’est pas stocké dans la contrainte : l’utilisateur l’écrit dans le NodeField de chargement, à la composante imposed (la duale propre de toute contrainte), au nœud-multiplicateur de la relation. Retrouver ce nœud à la main est fastidieux ; le helper le fait :

rhs = dirichlet.constraint_rhs([(noeud_contraint, u_d)])
rhs = mpc.constraint_rhs([(noeud_terme, g)])

• on désigne chaque relation par un nœud : le nœud contraint pour Dirichlet (un seul par relation), n’importe quel nœud-terme pour une MPC ;
• le helper résout ce nœud vers le nœud-multiplicateur de sa relation (via relations()) et y écrit la valeur, à la composante duale de la contrainte ;
• il renvoie un NodeField neuf sur tous les nœuds-multiplicateurs (les relations non citées valent 0), à fusionner dans le chargement global avec | : load | dirichlet.constraint_rhs(…) | mpc.constraint_rhs(…).

Le modèle passé doit porter exactement une contrainte (l’objet dirichlet ou mpc). Une erreur est levée si un nœud n’appartient à aucune relation, ou s’il en désigne plusieurs (ambigu).

Désigner par index de relation

Quand un même nœud participe à plusieurs relations (le keying par nœud est alors ambigu), on désigne la relation par son index (0-based, dans l’ordre de relations()) :

rhs = mpc.constraint_rhs_by_index([(index_relation, g)])

Le champ renvoyé et la fusion par | sont identiques ; une erreur est levée si un index dépasse le nombre de relations.

Contraintes disponibles

• Dirichlet — impose la valeur d’une variable primale (T = u_d, u_x = 0, …) sur un ensemble de nœuds. C’est la relation à un seul terme 1·u = u_d.
• Multi-points (MPC) — impose une relation linéaire à N termes Σₖ aₖ·u(nœudₖ, varₖ) = g entre plusieurs DOFs (égalité, périodicité, liaison affine…). Généralise Dirichlet.
• Baignage (embedded) — lie chaque nœud d’un maillage immergé à l’interpolation d’un maillage hôte au même point (u_c(p) = Σᵢ Nᵢ(ξ_p)·u_c(hôteᵢ)) : une barre baignée dans un volume. Les poids Nᵢ sont calculés par localisation de point à la construction ; c’est une MPC dont les coefficients varient par nœud.
• Contact (nœud-surface) — empêche les nœuds d’un maillage esclave de pénétrer une surface maître : une relation unilatérale (≥) par nœud esclave, à coefficients n·Nᵢ calculés par projection à la construction (petits glissements, sans frottement). Résolu par solve, qui y trouve des inégalités.

Une nouvelle contrainte suit le même patron : une struct implémentant SubModelKind, des blocs C/Cᵀ littéraux. Voir Ajouter une physique.