Calcul mécanique

Toujours la même plaque trouée : bord gauche encastré, effort réparti sur la moitié basse du trou (la « masse suspendue » du TD Cast3M). Trois volets, dans le même ordre que le support original — élasticité linéaire (sections 6-7-8), plasticité non linéaire pas à pas (section 9), contact unilatéral (section 10).

1. Élasticité linéaire

def construire_plaque_trouee():
    """A holed rectangular plate, built edge by edge with the
    mailleurs dédiés (`line`, `circle`), fusionnés en un seul
    submeshes by `pyrucast.mesh.consolidate` before `triangulate_surface` — as
    in `formation/maillage.py`. Also returns the submeshes the mechanics and
    the thermics need: left edge (clamped end), lower half of the hole
    (loading) and the whole hole (imposed temperature, reused as is to stay on
    the same nodes as `plaque`)."""
    coords = pc.Coords(2)
    p1 = coords.add_node([0.0, 0.0])
    p2 = coords.add_node([LONGUEUR, 0.0])
    p3 = coords.add_node([LONGUEUR, HAUTEUR])
    p4 = coords.add_node([0.0, HAUTEUR])

    bas = pc.mesh.line(p1, p2, 10)
    droit = pc.mesh.line(p2, p3, 4)
    haut = pc.mesh.line(p3, p4, 10)
    gauche = pc.mesh.line(p4, p1, 4)
    boucle_ext = pc.mesh.consolidate(bas | droit | haut | gauche)

    centre = coords.add_node(list(CENTRE_TROU))
    trou = pc.mesh.circle(centre, [0.0, 0.0, 1.0], RAYON_TROU, 16)

    # Outer loop CCW, hole clockwise (CW): the orientation
    # `triangulate_surface` expects (the hole is inverted, `trou` stays usable below).
    contour = boucle_ext | pc.mesh.invert(trou)
    plaque = pc.mesh.triangulate_surface(contour, "TRI3", size=0.02)

    # Lower half of the hole (y < centre): support of the hung mass's force,
    # like Cast3M's `PRES 'MASS'` over half the circle.
    y = pc.node_field.positions(trou, ["Y"])
    noeuds_bas_trou = pc.mesh.select(y, lt=CENTRE_TROU[1])
    arc_bas = pc.mesh.elements_on(trou, noeuds_bas_trou, strict=True)

    return coords, plaque, gauche, arc_bas, trou

Le modèle : élasticité (contraintes planes), encastrement u_x = u_y = 0 sur le bord gauche, effort réparti sur l’arc bas du trou (Cast3M FSUR 'MASS'/PRES 'MASS', ici pyrucast.model.flux en composante f_y) :

    encastrement = pc.mesh.to_poi1(gauche)
    multiplicateur = pc.mesh.translate(encastrement, [0.0, 0.0])

    modele = pc.model.elasticity(fes, "plane_stress")
    modele = modele | pc.model.dirichlet(modele, "u_x", encastrement, multiplicateur)
    modele = modele | pc.model.dirichlet(modele, "u_y", encastrement, multiplicateur)

    # The hung mass's force, spread over the lower half of the hole —
    # analogue de FSUR 'MASS' / PRES 'MASS' (Cast3M section 6).
    pression = -MASSE * G / (2.0 * 3.14159265 * RAYON_TROU)
    modele = modele | pc.model.flux(arc_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("alpha", ALPHA), ("phi_f_y", pression)]
    )
    effort = pc.node_field.external_forces(modele, materiaux)

    K = pc.matrix.stiffness(modele, materiaux)

    u1 = pc.solver.solve(K, effort)
    print(f"1) élasticité seule       : u_y(trou) ≈ {u1.min('u_y'):.6e} m")

2. + Dilatation thermique

On réutilise le champ de température résolu comme en Calcul thermique (température imposée sur le bord du trou) pour calculer une déformation thermique et l’ajouter au chargement — l’équivalent Cast3M EPTH :

    temperature = resoudre_thermique(plaque, trou)
    t_gauss = pc.element_field.interp_to_gauss(
        pc.node_field.restrict(temperature, plaque), fes
    )
    eps_th = pc.element_field.thermal_strain(t_gauss, materiaux, fes, T_REF)
    sig_th = pc.element_field.integrate_behavior(modele, eps_th, materiaux)
    f_th = (
        pc.node_field.divergence(sig_th, "sigma")
        .rename_component("div_sigma_x", "f_x")
        .rename_component("div_sigma_y", "f_y")
    )

    second_membre = f_th + pc.node_field.restrict_like(effort, f_th)
    u2 = pc.solver.solve(K, second_membre)
    print(f"2) + dilatation thermique : u_y(trou) ≈ {u2.min('u_y'):.6e} m")

Trois briques, à la main, sans opérateur « tout-en-un » — comme en Cast3M (EPTH + ELAS + BSIG) :

1. pyrucast.element_field.thermal_strain : \( \varepsilon_{\text{th}} = \alpha \cdot (T - T_{\text{ref}}) \), la même formule que Cast3M EPTH, à partir d’un champ de température aux points de Gauss (pyrucast.element_field.interp_to_gauss) ;
2. pyrucast.element_field.integrate_behavior : la pseudo-contrainte thermique \( \sigma_{\text{th}} = D : \varepsilon_{\text{th}} \) ;
3. pyrucast.node_field.internal_forces : la charge nodale équivalente \( F_{\text{th}} = \int B^T \sigma_{\text{th}} \, dV \) (Cast3M BSIG).

Le second membre se combine par addition de champs (+), pas par union (|) : internal_forces couvre tous les nœuds mécaniques, l’effort extérieur n’en couvre qu’une partie — pyrucast.node_field.restrict_like étend l’un sur le support de l’autre avant de les additionner. C’est la même distinction que la note de la page Calcul thermique : | pour des supports disjoints, +/- pour une véritable superposition sur un support commun.

Pour aller plus loin — matériau hétérogène (Cast3M section 8). Cast3M fait varier α(x) par une formule évaluée sur un champ aux points de Gauss (loi normale centrée sur la plaque). pyrucast le permettrait de la même façon : un ElementField accepte des valeurs non uniformes par (cellule, point de Gauss), et pyrucast.field.exp/l’arithmétique de champs (+ - * **) suffiraient à coder la formule — mais ce script ne le met pas en œuvre, faute d’un exemple testé à ce jour dans cette formation.

3. Plasticité parfaite — pas à pas

Le chargement dépasse maintenant la limite élastique. Comme Cast3M pilote ce cas par la procédure PASAPAS, pyrucast fournit pyrucast.thermomechanics.step_by_step : la boucle sur les pas de charge, un Newton modifié (rigidité élastique, réassemblée une fois par pas) et son accélération d’Anderson. Il suffit de remplacer model.elasticity par model.plasticity_perfect :

    encastrement = pc.mesh.to_poi1(bord_gauche)
    multiplicateur = pc.mesh.translate(encastrement, [0.0, 0.0])

    modele = pc.model.plasticity_perfect(fes, "plane_stress")
    modele = modele | pc.model.dirichlet(modele, "u_x", encastrement, multiplicateur)
    modele = modele | pc.model.dirichlet(modele, "u_y", encastrement, multiplicateur)

L’historique de charge est un Evolution à valeur champ, interpolée linéairement en pseudo-temps — Cast3M EVOL 'MANU' :

    pression = -FACTEUR_CHARGE * MASSE * G / (2.0 * 3.14159265 * RAYON_TROU)
    modele = modele | pc.model.flux(arc_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("sigma_y", SIGMA_Y), ("phi_f_y", pression)]
    )
    effort_final = pc.node_field.external_forces(modele, materiaux)
    charge = pc.Evolution(
        [(0.0, effort_final * 0.0), (1.0, effort_final)], out_of_range="clamp"
    )

    # Free DOFs (outside the clamped end) to norm the Newton residual —
    # without which the large support reactions mask the real convergence.
    x = pc.node_field.positions(plaque, ["X"])
    ddl_libres = pc.mesh.select(x, gt=1e-6)

    data = {
        "times": [0.0, 0.2, 0.4, 0.55, 0.7],  # pseudo-temps ∈ [0, 1]
        "model": modele,
        "loads": charge,
        "materials": materiaux,
        "free_mesh": ddl_libres,
        "max_newton": 200,
    }
    pc.thermomechanics.step_by_step(data)

Piège pyrucast. Sans free_mesh, la norme du résidu de Newton porte sur tous les nœuds, y compris les nœuds encastrés — dont la réaction d’appui, potentiellement énorme, empêche toute convergence. free_mesh restreint la norme aux degrés de liberté réellement libres (ici : tous les nœuds d’abscisse x > 0). C’est l’équivalent, en plus explicite, du traitement automatique des blocages par Cast3M dans RESO/PASAPAS.

Non disponible dans pyrucast. model.plasticity_perfect ne consomme pas encore la composante matériau optionnelle alpha — la dépendance de sigma_y à la température (Cast3M section 9.2 : EVOL 'MANU' 'T' ... 'SIGY' ...) n’a donc pas d’équivalent testé ici ; seule la plasticité isotherme est couverte.

4. Contact unilatéral

Cast3M pilote aussi le contact par PASAPAS (table tab3, section 10). pyrucast ne compose pas encore thermique + plasticité + contact dans un même appel step_by_step ; le contact se résout directement par le solveur actif-set pyrucast.solver.solve_unilateral, sur un patch-test classique (deux blocs superposés, jeu initial, pression sur le bloc du haut) :

    coords = pc.Coords(2)
    mesh_bas, bas = bloc(coords, 0.0)
    mesh_haut, haut = bloc(coords, 1.0 + G0)
    mesh = mesh_bas | mesh_haut
    fes = pc.FiniteElementSpace(mesh)

    # Master: top edge of the lower block (`contour` already orients the
    # boundary counter-clockwise, so this edge naturally runs right to left —
    # the associated normal points towards +y). Slave: nodes of the upper
    # block's bottom edge.
    maitre = bord_horizontal(mesh_bas, 1.0)
    esclave = pc.mesh.poi1_from_nodes([haut[idx(i, 0)] for i in range(N + 1)])

    elasticite = pc.model.elasticity(fes, "plane_stress")
    contact = pc.model.contact(elasticite, esclave, maitre, ["u_x", "u_y"])

    modele = pc.model.elasticity(fes, "plane_stress")
    modele = modele | clamp(modele, bas + haut, "u_x")
    modele = modele | clamp(modele, [bas[idx(i, 0)] for i in range(N + 1)], "u_y")
    modele = modele | contact

    bord_haut = bord_horizontal(mesh_haut, 2.0 + G0)
    bord_haut_fes = pc.FiniteElementSpace(bord_haut)
    modele = modele | pc.model.flux(bord_haut_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", 0.0), ("phi_f_y", -S)]
    )
    traction = pc.node_field.external_forces(modele, materiaux)

    # `contact_gaps()` supplies the contact's right-hand side — the Cast3M
    # equivalent of preparing the unilateral problem before RESO.
    second_membre = traction | modele.contact_gaps()

    K = pc.matrix.stiffness(modele, materiaux)
    solution = pc.solver.solve_unilateral(K, modele, second_membre)

model.contact_gaps() fournit le second membre associé au contact — jeu initial à combler avant que les multiplicateurs lambda_contact ne portent une réaction non nulle. Voir Contact (nœud-surface) pour le détail mathématique (formulation active-set, jeux, multiplicateurs).

Scripts complets

"""Formation débutant — 3. Calcul mécanique (élasticité linéaire).

Reprend la plaque trouée : bord gauche encastré, effort ponctuel (masse
mass) spread over the lower half of the hole. Three load cases follow one
another, like sections 6/7/8 of the Cast3M training:

1. **élasticité pure** — effort seul ;
2. **+ dilatation thermique** — on réutilise le champ de température de
   `formation/thermique.py` (`ε_th = α·(T − T_ref)`, opérateur
   `field.thermal_strain`, l'équivalent Cast3M `EPTH`) ;
3. a paragraph (no code tested here) on the **heterogeneous material**
   (Cast3M varies `alpha(x)` through a formula on a field at the Gauss
   points) — see the book page for the detail.

Lancement ::

    maturin develop --release
    python formation/mecanique.py

    # To regenerate the book figure (book/src/formation/img/):
    # PYRUCAST_FORMATION_IMG_DIR=book/src/formation/img python formation/mecanique.py
"""

import os
import tempfile

import pyrucast as pc

LONGUEUR, HAUTEUR = 0.30, 0.10  # m
RAYON_TROU = 0.025  # m
CENTRE_TROU = (0.75 * LONGUEUR, HAUTEUR / 2.0)

E, NU, ALPHA = 200e9, 0.3, 1e-5  # acier
MASSE, G = 2500.0, 9.81  # kg, m/s^2 — mass hung from the hole
T_REF, T_IMPOSEE = 20.0, 250.0  # °C — dilatation thermique
K_COND = 50.0  # W/m/K


# ANCHOR: construction
def construire_plaque_trouee():
    """A holed rectangular plate, built edge by edge with the
    mailleurs dédiés (`line`, `circle`), fusionnés en un seul
    submeshes by `pyrucast.mesh.consolidate` before `triangulate_surface` — as
    in `formation/maillage.py`. Also returns the submeshes the mechanics and
    the thermics need: left edge (clamped end), lower half of the hole
    (loading) and the whole hole (imposed temperature, reused as is to stay on
    the same nodes as `plaque`)."""
    coords = pc.Coords(2)
    p1 = coords.add_node([0.0, 0.0])
    p2 = coords.add_node([LONGUEUR, 0.0])
    p3 = coords.add_node([LONGUEUR, HAUTEUR])
    p4 = coords.add_node([0.0, HAUTEUR])

    bas = pc.mesh.line(p1, p2, 10)
    droit = pc.mesh.line(p2, p3, 4)
    haut = pc.mesh.line(p3, p4, 10)
    gauche = pc.mesh.line(p4, p1, 4)
    boucle_ext = pc.mesh.consolidate(bas | droit | haut | gauche)

    centre = coords.add_node(list(CENTRE_TROU))
    trou = pc.mesh.circle(centre, [0.0, 0.0, 1.0], RAYON_TROU, 16)

    # Outer loop CCW, hole clockwise (CW): the orientation
    # `triangulate_surface` expects (the hole is inverted, `trou` stays usable below).
    contour = boucle_ext | pc.mesh.invert(trou)
    plaque = pc.mesh.triangulate_surface(contour, "TRI3", size=0.02)

    # Lower half of the hole (y < centre): support of the hung mass's force,
    # like Cast3M's `PRES 'MASS'` over half the circle.
    y = pc.node_field.positions(trou, ["Y"])
    noeuds_bas_trou = pc.mesh.select(y, lt=CENTRE_TROU[1])
    arc_bas = pc.mesh.elements_on(trou, noeuds_bas_trou, strict=True)

    return coords, plaque, gauche, arc_bas, trou


# ANCHOR_END: construction


def resoudre_thermique(plaque, trou):
    """Ré-sout la thermique de `formation/thermique.py` (version simplifiée,
    without convection or source, just T imposed on the hole) so as to reuse
    a non-uniform field in the second load case below.

    Important: we reuse the `trou` returned by `construire_plaque_trouee` —
    hence the same nodes as the hole's edge in `plaque` — rather than
    rebuilding a separate circle, which would give nodes disjoint from the
    real mesh and a Dirichlet with no effect on the solution."""
    fes = pc.FiniteElementSpace(plaque)
    modele_th = pc.model.heat_conduction(fes)

    trou_poi1 = pc.mesh.to_poi1(trou)
    multiplicateur = pc.mesh.translate(trou_poi1, [0.0, 0.0])
    modele_th = modele_th | pc.model.dirichlet(
        modele_th, "T", trou_poi1, multiplicateur
    )

    materiaux_th = pc.element_field.material_field(modele_th, [("k", K_COND)])
    temperature_imposee = pc.NodeField(multiplicateur, ["imposed_T"])
    temperature_imposee[0].add_to_component("imposed_T", T_IMPOSEE)

    K_th = pc.matrix.stiffness(modele_th, materiaux_th)
    return pc.solver.solve(K_th, temperature_imposee)


def main() -> None:
    _coords, plaque, gauche, arc_bas, trou = construire_plaque_trouee()
    fes = pc.FiniteElementSpace(plaque)
    arc_fes = pc.FiniteElementSpace(arc_bas)

    # ANCHOR: modele_elastique
    encastrement = pc.mesh.to_poi1(gauche)
    multiplicateur = pc.mesh.translate(encastrement, [0.0, 0.0])

    modele = pc.model.elasticity(fes, "plane_stress")
    modele = modele | pc.model.dirichlet(modele, "u_x", encastrement, multiplicateur)
    modele = modele | pc.model.dirichlet(modele, "u_y", encastrement, multiplicateur)

    # The hung mass's force, spread over the lower half of the hole —
    # analogue de FSUR 'MASS' / PRES 'MASS' (Cast3M section 6).
    pression = -MASSE * G / (2.0 * 3.14159265 * RAYON_TROU)
    modele = modele | pc.model.flux(arc_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("alpha", ALPHA), ("phi_f_y", pression)]
    )
    effort = pc.node_field.external_forces(modele, materiaux)

    K = pc.matrix.stiffness(modele, materiaux)
    # ANCHOR_END: modele_elastique

    # ANCHOR: cas1_elastique
    u1 = pc.solver.solve(K, effort)
    print(f"1) élasticité seule       : u_y(trou) ≈ {u1.min('u_y'):.6e} m")
    # ANCHOR_END: cas1_elastique

    # ANCHOR: cas2_thermique
    temperature = resoudre_thermique(plaque, trou)
    t_gauss = pc.element_field.interp_to_gauss(
        pc.node_field.restrict(temperature, plaque), fes
    )
    eps_th = pc.element_field.thermal_strain(t_gauss, materiaux, fes, T_REF)
    sig_th = pc.element_field.integrate_behavior(modele, eps_th, materiaux)
    f_th = (
        pc.node_field.divergence(sig_th, "sigma")
        .rename_component("div_sigma_x", "f_x")
        .rename_component("div_sigma_y", "f_y")
    )

    second_membre = f_th + pc.node_field.restrict_like(effort, f_th)
    u2 = pc.solver.solve(K, second_membre)
    print(f"2) + dilatation thermique : u_y(trou) ≈ {u2.min('u_y'):.6e} m")
    # ANCHOR_END: cas2_thermique

    # u2 also carries the Dirichlet's Lagrange multipliers: only (u_x, u_y)
    # are kept before computing a strain.
    u2_propre = pc.node_field.restrict_like(u2, pc.NodeField(plaque, ["u_x", "u_y"]))
    contraintes = pc.element_field.integrate_behavior(
        modele, pc.element_field.deformation(u2_propre, fes) - eps_th, materiaux
    )
    print(f"   σ_xx max ≈ {contraintes.max('sigma_xx'):.3e} Pa")

    out = os.environ.get("PYRUCAST_FORMATION_IMG_DIR", tempfile.gettempdir())
    chemin = os.path.join(out, "mecanique-deplacement.svg")
    plaque.plot(save=chemin, field=u2, component="u_y", cmap="coolwarm", smooth=1)
    print(f"Displacement u_y written to {chemin}")


if __name__ == "__main__":
    main()

"""Formation débutant — 4. Mécanique non linéaire (plasticité).

Picks the holed plate clamped on the left back up, with the hung mass's force
**ramped up** until it passes the elastic limit —
l'équivalent Python de la table `PASAPAS` de Cast3M (section 9 de la
formation) : ``pyrucast.thermomechanics.step_by_step`` orchestre la boucle
over the load steps and, at each step, a **modified** Newton (elastic
stiffness, sped up by Anderson acceleration).

Replacing ``model.elasticity`` with ``Model.plasticity`` in
`formation/mecanique.py` is all it takes to get this script — the same call
``step_by_step`` gère la boucle non linéaire.

A pyrucast caveat, specific to this version of the library: plasticity
(`Model.plasticity`) does not consume the
matériau optionnelle `alpha` (dilatation thermique, Cast3M `EPTH`) — le
thermo-plastic coupling of Cast3M's section 9.2 (where `sigma_y` depends on
temperature) is therefore not covered here.

Lancement ::

    maturin develop --release
    python formation/plasticite.py

    # To regenerate the book figure (book/src/formation/img/):
    # PYRUCAST_FORMATION_IMG_DIR=book/src/formation/img python formation/plasticite.py
"""

import os
import tempfile

import pyrucast as pc

LONGUEUR, HAUTEUR = 0.30, 0.10  # m
RAYON_TROU = 0.025  # m
CENTRE_TROU = (0.75 * LONGUEUR, HAUTEUR / 2.0)

E, NU = 200e9, 0.3
# σy deliberately modest: this training's geometry and loading are not at the
# scale of a real steel — what matters is to bring out a plastic zone in a
# few steps, not the material's physical reality.

SIGMA_Y = 5e6
MASSE, G = 2500.0, 9.81
FACTEUR_CHARGE = 6.0  # multiplier on the hung mass, to pass σy


def construire_plaque_trouee():
    """Same geometry as `formation/mecanique.py`."""
    coords = pc.Coords(2)
    p1 = coords.add_node([0.0, 0.0])
    p2 = coords.add_node([LONGUEUR, 0.0])
    p3 = coords.add_node([LONGUEUR, HAUTEUR])
    p4 = coords.add_node([0.0, HAUTEUR])

    bas = pc.mesh.line(p1, p2, 10)
    droit = pc.mesh.line(p2, p3, 4)
    haut = pc.mesh.line(p3, p4, 10)
    bord_gauche = pc.mesh.line(p4, p1, 4)
    boucle_ext = pc.mesh.consolidate(bas | droit | haut | bord_gauche)

    centre = coords.add_node(list(CENTRE_TROU))
    trou = pc.mesh.circle(centre, [0.0, 0.0, 1.0], RAYON_TROU, 16)

    # Outer loop CCW, hole clockwise (CW): the orientation
    # `triangulate_surface` expects (the hole is inverted, `trou` stays usable below).
    contour = boucle_ext | pc.mesh.invert(trou)
    plaque = pc.mesh.triangulate_surface(contour, "TRI3", size=0.02)

    y = pc.node_field.positions(trou, ["Y"])
    noeuds_bas_trou = pc.mesh.select(y, lt=CENTRE_TROU[1])
    arc_bas = pc.mesh.elements_on(trou, noeuds_bas_trou, strict=True)

    return coords, plaque, bord_gauche, arc_bas


def main() -> None:
    _coords, plaque, bord_gauche, arc_bas = construire_plaque_trouee()
    fes = pc.FiniteElementSpace(plaque)
    arc_fes = pc.FiniteElementSpace(arc_bas)

    # ANCHOR: modele_plastique
    encastrement = pc.mesh.to_poi1(bord_gauche)
    multiplicateur = pc.mesh.translate(encastrement, [0.0, 0.0])

    modele = pc.model.plasticity_perfect(fes, "plane_stress")
    modele = modele | pc.model.dirichlet(modele, "u_x", encastrement, multiplicateur)
    modele = modele | pc.model.dirichlet(modele, "u_y", encastrement, multiplicateur)
    # ANCHOR_END: modele_plastique

    # ANCHOR: chargement_evolution
    pression = -FACTEUR_CHARGE * MASSE * G / (2.0 * 3.14159265 * RAYON_TROU)
    modele = modele | pc.model.flux(arc_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", NU), ("sigma_y", SIGMA_Y), ("phi_f_y", pression)]
    )
    effort_final = pc.node_field.external_forces(modele, materiaux)
    charge = pc.Evolution(
        [(0.0, effort_final * 0.0), (1.0, effort_final)], out_of_range="clamp"
    )
    # ANCHOR_END: chargement_evolution

    # ANCHOR: pas_a_pas
    # Free DOFs (outside the clamped end) to norm the Newton residual —
    # without which the large support reactions mask the real convergence.
    x = pc.node_field.positions(plaque, ["X"])
    ddl_libres = pc.mesh.select(x, gt=1e-6)

    data = {
        "times": [0.0, 0.2, 0.4, 0.55, 0.7],  # pseudo-temps ∈ [0, 1]
        "model": modele,
        "loads": charge,
        "materials": materiaux,
        "free_mesh": ddl_libres,
        "max_newton": 200,
    }
    pc.thermomechanics.step_by_step(data)
    # ANCHOR_END: pas_a_pas

    print(f"{'t':>6} {'itérations':>11} {'anderson':>9} {'convergé':>9} {'p_max':>12}")
    for r in data["results"]:
        p_max = r["state"].max("p") if r["state"] is not None else 0.0
        print(
            f"{r['time']:>6.2f} {r['mech_iters']:>11} {r['mech_anderson']:>9} "
            f"{r['converged']!s:>9} {p_max:>12.3e}"
        )

    dernier = data["results"][-1]
    print(f"\nzone plastique développée : p_max = {dernier['state'].max('p'):.3e}")

    out = os.environ.get("PYRUCAST_FORMATION_IMG_DIR", tempfile.gettempdir())
    chemin = os.path.join(out, "plasticite.svg")
    plaque.plot(
        save=chemin, field=dernier["state"], component="p", cmap="viridis", smooth=0
    )
    print(f"Plastic zone (p) written to {chemin}")


if __name__ == "__main__":
    main()

"""Formation débutant — 5. Contact (unilatéral, nœud-surface).

A classic patch test: two elastic blocks stacked along `y`, separated by an
initial gap `G0`. A pressure on the upper block closes the contact and
transmits a uniform stress across the interface — the pyrucast equivalent of
Cast3M's node-to-surface contact (section 10 of the training), driven here
straight by the active-set solver `solve_unilateral` rather than by
`step_by_step` (which cannot yet compose thermics, plasticity and contact in
a single table).

Lancement ::

    maturin develop --release
    python formation/contact.py

    # To regenerate the book figure (book/src/formation/img/):
    # PYRUCAST_FORMATION_IMG_DIR=book/src/formation/img python formation/contact.py
"""

import os
import tempfile

import pyrucast as pc

E = 100.0
S = 5.0  # pression appliquée
G0 = 0.01  # initial gap between the two blocks
N = 2  # N×N grid of QUA4 per block


def idx(i, j):
    return j * (N + 1) + i


def bloc(coords: pc.Coords, y0: float):
    """Bloc `[0,1] × [y0, y0+1]`, grille N×N de QUA4 — mailleurs dédiés
    (`line` for the bottom/top edges, `sweep` between them, as in
    `formation/maillage.py`). Returns `(mesh, grille)`, `grille[idx(i,j)]`
    étant le nœud `(i,j)` (`i` : abscisse, `j` : ordonnée)."""
    bas = pc.mesh.line(coords.add_node([0.0, y0]), coords.add_node([1.0, y0]), N)
    haut = pc.mesh.line(
        coords.add_node([0.0, y0 + 1.0]), coords.add_node([1.0, y0 + 1.0]), N
    )
    mesh = pc.mesh.sweep(bas, haut, N)

    grille = [None] * ((N + 1) * (N + 1))
    for cy in range(N):
        for cx in range(N):
            cell = cy * N + cx
            grille[idx(cx, cy)] = mesh.node(0, cell, 0)
            grille[idx(cx + 1, cy)] = mesh.node(0, cell, 1)
            grille[idx(cx + 1, cy + 1)] = mesh.node(0, cell, 2)
            grille[idx(cx, cy + 1)] = mesh.node(0, cell, 3)
    return mesh, grille


def clamp(target, nodes, var):
    imposed = pc.mesh.poi1_from_nodes(nodes)
    multiplier = pc.mesh.barycenter(imposed)
    return pc.model.dirichlet(target, var, imposed, multiplier)


def bord_horizontal(mesh: pc.Mesh, y: float) -> pc.Mesh:
    """Extracts, among `mesh`'s border segments (`pyrucast.mesh.border`, the
    Cast3M `CONTOUR` equivalent), those at ordinate `y` — an existing edge of
    the mesh, not a line rebuilt beside it (`line` would make
    nouveaux nœuds, disjoints de `mesh`)."""
    frontiere = pc.mesh.border(mesh)
    ordonnee = pc.node_field.positions(frontiere, ["Y"])
    noeuds = pc.mesh.select(ordonnee, ge=y - 1e-9, le=y + 1e-9)
    return pc.mesh.elements_on(frontiere, noeuds, strict=True)


def main() -> None:
    # ANCHOR: geometrie_contact
    coords = pc.Coords(2)
    mesh_bas, bas = bloc(coords, 0.0)
    mesh_haut, haut = bloc(coords, 1.0 + G0)
    mesh = mesh_bas | mesh_haut
    fes = pc.FiniteElementSpace(mesh)

    # Master: top edge of the lower block (`contour` already orients the
    # boundary counter-clockwise, so this edge naturally runs right to left —
    # the associated normal points towards +y). Slave: nodes of the upper
    # block's bottom edge.
    maitre = bord_horizontal(mesh_bas, 1.0)
    esclave = pc.mesh.poi1_from_nodes([haut[idx(i, 0)] for i in range(N + 1)])

    elasticite = pc.model.elasticity(fes, "plane_stress")
    contact = pc.model.contact(elasticite, esclave, maitre, ["u_x", "u_y"])
    # ANCHOR_END: geometrie_contact

    # ANCHOR: modele_contact
    modele = pc.model.elasticity(fes, "plane_stress")
    modele = modele | clamp(modele, bas + haut, "u_x")
    modele = modele | clamp(modele, [bas[idx(i, 0)] for i in range(N + 1)], "u_y")
    modele = modele | contact

    # ANCHOR_END: modele_contact

    # ANCHOR: chargement_contact
    bord_haut = bord_horizontal(mesh_haut, 2.0 + G0)
    bord_haut_fes = pc.FiniteElementSpace(bord_haut)
    modele = modele | pc.model.flux(bord_haut_fes, modele, "f_y")
    materiaux = pc.element_field.material_field(
        modele, [("E", E), ("nu", 0.0), ("phi_f_y", -S)]
    )
    traction = pc.node_field.external_forces(modele, materiaux)

    # `contact_gaps()` supplies the contact's right-hand side — the Cast3M
    # equivalent of preparing the unilateral problem before RESO.
    second_membre = traction | modele.contact_gaps()
    # ANCHOR_END: chargement_contact

    # ANCHOR: resolution_contact
    K = pc.matrix.stiffness(modele, materiaux)
    solution = pc.solver.solve_unilateral(K, modele, second_membre)
    # ANCHOR_END: resolution_contact

    print(f"Pression appliquée : {S}")
    for j in range(N + 1):
        uy_bas = solution.value(bas[idx(0, j)], "u_y")
        uy_haut = solution.value(haut[idx(0, j)], "u_y")
        print(f"  y={j / N:.2f} : u_y(bas)={uy_bas:.6e}  u_y(haut)={uy_haut:.6e}")

    # Réactions de contact : Σ(−λᵢ) doit reconstituer l'effort appliqué S.
    maillage_mult = contact.multiplier_mesh()
    lambdas = [
        solution.value(maillage_mult.node(0, r, 0), "lambda_contact")
        for r in range(N + 1)
    ]
    print(f"\nΣ(−λ) = {sum(-lam for lam in lambdas):.6f}  (attendu {S})")

    out = os.environ.get("PYRUCAST_FORMATION_IMG_DIR", tempfile.gettempdir())
    chemin = os.path.join(out, "contact.svg")
    maillage_mult.plot(
        save=chemin, field=solution, component="lambda_contact", cmap="viridis"
    )
    print(f"Contact reaction (λ) written to {chemin}")


if __name__ == "__main__":
    main()

Suite : Compléments — éléments structuraux et export des résultats.