19 — Elastoplastic Pushover (Uniaxial Bar)¶

Curriculum slot: Tier 5, slot 19.
Prerequisite: 06 — Sections Catalog, 18 — Buckling.
Strategy used for results: spec.capture(...) — multi-step nonlinear with broad coverage.

Purpose¶

Pushover is displacement-controlled nonlinear static analysis — the standard tool for capacity curves (base reaction vs controlled displacement) up to / beyond yield. Two ingredients new to the curriculum:

1. A nonlinear material law. We use Steel01 (elastic-perfectly-plastic, no hardening) so the analytical reference is bilinear and clean.
2. Displacement control as the integrator. DisplacementControl ramps a target DOF and the analysis returns the reaction.

Problem¶

A 1-D bar of length $L$ and cross-section $A$ — Steel01 with $E$, $f_y$, $b = 0$. Left end fixed; right end pulled by displacement $u$. The reaction at the fixed end traces:

$$ F(u) \;=\; \begin{cases} E\,A\,u/L, & u < u_y \\ f_y\,A = F_y, & u \geq u_y \end{cases}, \quad u_y = \dfrac{f_y\,L}{E}. $$

Why this notebook is interesting¶

• The first multi-step nonlinear analysis of the curriculum — the recorder fires at every load step instead of once at the end.
• spec.capture shines: time-history of displacement_x at the pulled node + reaction_force_x at the fixed end gives the capacity curve without any manual nodeDisp / nodeReaction loops in the analysis driver.
• The plot is a single-step expression on the slabs: plt.plot(u_slab.values, F_slab.values).

1. Imports and parameters¶

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import openseespy.opensees as ops

from apeGmsh import apeGmsh, Results
from apeGmsh.results.spec import Recorders

# Geometry / material (consistent units: m, Pa, N)
L  = 1.0
A  = 1.0e-4
E  = 2.0e11
fy = 300.0e6
Fy = fy * A         # force at yield
uy = fy * L / E      # displacement at yield

# Displacement-controlled ramp 0 → 4 u_y
N_STEPS = 80
u_max   = 4.0 * uy
dU      = u_max / N_STEPS

# Mesh
N_ELEM = 10
LC = L / N_ELEM

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import openseespy.opensees as ops from apeGmsh import apeGmsh, Results from apeGmsh.results.spec import Recorders # Geometry / material (consistent units: m, Pa, N) L = 1.0 A = 1.0e-4 E = 2.0e11 fy = 300.0e6 Fy = fy * A # force at yield uy = fy * L / E # displacement at yield # Displacement-controlled ramp 0 → 4 u_y N_STEPS = 80 u_max = 4.0 * uy dU = u_max / N_STEPS # Mesh N_ELEM = 10 LC = L / N_ELEM

2. Geometry + mesh¶

In [ ]:

g_ctx = apeGmsh(model_name="19_pushover", verbose=False)
g = g_ctx.__enter__()

p_L = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC)
p_R = g.model.geometry.add_point(L,   0.0, 0.0, lc=LC)
ln  = g.model.geometry.add_line(p_L, p_R)
g.model.sync()

g.physical.add(0, [p_L], name="left")
g.physical.add(0, [p_R], name="right")
g.physical.add(1, [ln],  name="bar")

g.mesh.structured.set_transfinite_curve(ln, N_ELEM + 1)
g.mesh.generation.generate(1)
fem = g.mesh.queries.get_fem_data()
print(f"mesh: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

g_ctx = apeGmsh(model_name="19_pushover", verbose=False) g = g_ctx.__enter__() p_L = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC) p_R = g.model.geometry.add_point(L, 0.0, 0.0, lc=LC) ln = g.model.geometry.add_line(p_L, p_R) g.model.sync() g.physical.add(0, [p_L], name="left") g.physical.add(0, [p_R], name="right") g.physical.add(1, [ln], name="bar") g.mesh.structured.set_transfinite_curve(ln, N_ELEM + 1) g.mesh.generation.generate(1) fem = g.mesh.queries.get_fem_data() print(f"mesh: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

3. Build the OpenSees model — vanilla openseespy¶

1-D model (-ndm 1 -ndf 1), Steel01 material with b = 0 (perfectly plastic), and truss elements. Left end fixed; the right end will be pulled by DisplacementControl. The reference load (a unit pull) defines the controlled DOF — its magnitude doesn't matter.

In [ ]:

ops.wipe()
ops.model("basic", "-ndm", 1, "-ndf", 1)

for nid, xyz in fem.nodes.get():
    ops.node(int(nid), float(xyz[0]))

MAT_TAG = 1
ops.uniaxialMaterial("Steel01", MAT_TAG, fy, E, 0.0)

for group in fem.elements.get(target="bar"):
    for eid, conn in zip(group.ids, group.connectivity):
        ops.element(
            "truss", int(eid),
            int(conn[0]), int(conn[1]),
            A, MAT_TAG,
        )

left_id  = int(next(iter(fem.nodes.get(target="left").ids)))
right_id = int(next(iter(fem.nodes.get(target="right").ids)))
ops.fix(left_id, 1)

# Reference load — defines the direction; magnitude irrelevant
ops.timeSeries("Linear", 1)
ops.pattern("Plain", 1, 1)
ops.load(right_id, 1.0)

ops.wipe() ops.model("basic", "-ndm", 1, "-ndf", 1) for nid, xyz in fem.nodes.get(): ops.node(int(nid), float(xyz[0])) MAT_TAG = 1 ops.uniaxialMaterial("Steel01", MAT_TAG, fy, E, 0.0) for group in fem.elements.get(target="bar"): for eid, conn in zip(group.ids, group.connectivity): ops.element( "truss", int(eid), int(conn[0]), int(conn[1]), A, MAT_TAG, ) left_id = int(next(iter(fem.nodes.get(target="left").ids))) right_id = int(next(iter(fem.nodes.get(target="right").ids))) ops.fix(left_id, 1) # Reference load — defines the direction; magnitude irrelevant ops.timeSeries("Linear", 1) ops.pattern("Plain", 1, 1) ops.load(right_id, 1.0)

4. Declare recorders¶

• displacement_x at the controlled (right) end.
• reaction_force_x at the fixed (left) end.

Both record at every step — the recorder fires once per successful analyze(1), so a single capture stage gives us the full capacity curve.

In [ ]:

recorders = Recorders()

recorders.nodes(
    components=["displacement_x"], pg="right", name="u_right",
)
recorders.nodes(
    components=["reaction_force_x"], pg="left", name="R_left",
)

spec = recorders.resolve(fem, ndm=1, ndf=1)
print(spec)

recorders = Recorders() recorders.nodes( components=["displacement_x"], pg="right", name="u_right", ) recorders.nodes( components=["reaction_force_x"], pg="left", name="R_left", ) spec = recorders.resolve(fem, ndm=1, ndf=1) print(spec)

5. Run the pushover with spec.capture(...)¶

Standard nonlinear recipe — Newton algorithm, displacement control with step dU. After each successful analyze(1) we call cap.step(t=...) so the capture writes one chunk per load step. Setting t = (k+1) * dU makes the time axis equal to the target displacement — convenient when reading back.

In [ ]:

from apeGmsh import workdir
OUT = workdir()
results_path = OUT / "capture.h5"
if results_path.exists():
    results_path.unlink()

with spec.capture(results_path, fem=fem, ndm=1, ndf=1) as cap:
    cap.begin_stage("pushover", kind="static")

    ops.system("BandGeneral")
    ops.numberer("Plain")
    ops.constraints("Plain")
    ops.test("NormDispIncr", 1e-10, 20)
    ops.algorithm("Newton")
    ops.integrator("DisplacementControl", right_id, 1, dU)
    ops.analysis("Static")

    converged = 0
    for k in range(N_STEPS):
        status = ops.analyze(1)
        if status != 0:
            print(f"  diverged at step {k + 1}")
            break
        cap.step(t=(k + 1) * dU)
        converged += 1

    cap.end_stage()

print(f"converged {converged} / {N_STEPS} steps")
print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

from apeGmsh import workdir OUT = workdir() results_path = OUT / "capture.h5" if results_path.exists(): results_path.unlink() with spec.capture(results_path, fem=fem, ndm=1, ndf=1) as cap: cap.begin_stage("pushover", kind="static") ops.system("BandGeneral") ops.numberer("Plain") ops.constraints("Plain") ops.test("NormDispIncr", 1e-10, 20) ops.algorithm("Newton") ops.integrator("DisplacementControl", right_id, 1, dU) ops.analysis("Static") converged = 0 for k in range(N_STEPS): status = ops.analyze(1) if status != 0: print(f" diverged at step {k + 1}") break cap.step(t=(k + 1) * dU) converged += 1 cap.end_stage() print(f"converged {converged} / {N_STEPS} steps") print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

6. Read results back — capacity curve from slabs¶

Two queries — pulled-end displacement and fixed-end reaction — and the capacity curve falls out as numpy arrays without any manual nodeDisp / nodeReaction loop in user code.

In [ ]:

results = Results.from_native(results_path, fem=fem)
print(results.inspect.summary())

results = Results.from_native(results_path, fem=fem) print(results.inspect.summary())

In [ ]:

u_slab = results.nodes.get(component="displacement_x", pg="right")
R_slab = results.nodes.get(component="reaction_force_x", pg="left")

u_hist = u_slab.values[:, 0]                # shape (T,)
F_hist = -R_slab.values.sum(axis=1)         # ON the bar at the left = − reaction
print(f"steps recorded: {len(u_hist)}")
print(f"u_max          : {u_hist[-1]:.4e}  m")
print(f"F at u_max     : {F_hist[-1]:.4e}  N")

u_slab = results.nodes.get(component="displacement_x", pg="right") R_slab = results.nodes.get(component="reaction_force_x", pg="left") u_hist = u_slab.values[:, 0] # shape (T,) F_hist = -R_slab.values.sum(axis=1) # ON the bar at the left = − reaction print(f"steps recorded: {len(u_hist)}") print(f"u_max : {u_hist[-1]:.4e} m") print(f"F at u_max : {F_hist[-1]:.4e} N")

7. Verification¶

Three checks against the closed-form bilinear capacity curve:

• Pre-yield slope $\to E$ (linear regression on $u < 0.8 u_y$).
• Yield plateau — mean of last 20% of steps should equal $F_y$.
• Yield onset — first step where $F \ge 0.99 F_y$ should occur near $u = u_y$.

In [ ]:

elastic_mask = u_hist < 0.8 * uy
slope = np.polyfit(u_hist[elastic_mask], F_hist[elastic_mask], 1)[0]
E_recovered = slope * L / A

tail = F_hist[int(0.8 * len(F_hist)):]
F_plateau = float(np.mean(tail))

idx_yield  = int(next(i for i, f in enumerate(F_hist) if f >= 0.99 * Fy))
u_at_yield = u_hist[idx_yield]

print("Pre-yield modulus")
print(f"  recovered E :  {E_recovered:.4e}  Pa")
print(f"  reference E :  {E:.4e}  Pa")
print(f"  error       :  {abs(E_recovered - E)/E*100:.4f} %")
print()
print("Yield plateau")
print(f"  mean(tail)   :  {F_plateau:.4e}  N")
print(f"  Fy = fy A    :  {Fy:.4e}  N")
print(f"  error        :  {abs(F_plateau - Fy)/Fy*100:.4f} %")
print()
print("Yield onset")
print(f"  first step F ≥ 0.99 Fy at u = {u_at_yield:.4e}  m")
print(f"  uy = fy L / E              : {uy:.4e}  m")

elastic_mask = u_hist < 0.8 * uy slope = np.polyfit(u_hist[elastic_mask], F_hist[elastic_mask], 1)[0] E_recovered = slope * L / A tail = F_hist[int(0.8 * len(F_hist)):] F_plateau = float(np.mean(tail)) idx_yield = int(next(i for i, f in enumerate(F_hist) if f >= 0.99 * Fy)) u_at_yield = u_hist[idx_yield] print("Pre-yield modulus") print(f" recovered E : {E_recovered:.4e} Pa") print(f" reference E : {E:.4e} Pa") print(f" error : {abs(E_recovered - E)/E*100:.4f} %") print() print("Yield plateau") print(f" mean(tail) : {F_plateau:.4e} N") print(f" Fy = fy A : {Fy:.4e} N") print(f" error : {abs(F_plateau - Fy)/Fy*100:.4f} %") print() print("Yield onset") print(f" first step F ≥ 0.99 Fy at u = {u_at_yield:.4e} m") print(f" uy = fy L / E : {uy:.4e} m")

8. Plot the capacity curve¶

In [ ]:

u_fine = np.linspace(0, u_hist[-1], 200)
F_closed = np.where(u_fine < uy, E * A * u_fine / L, Fy)

fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(u_fine * 1e3, F_closed * 1e-3, "-", color="k", lw=1, label="closed-form bilinear")
ax.plot(u_hist * 1e3, F_hist * 1e-3, "o", ms=3, color="C0", label="FEM (capture)")
ax.axvline(uy * 1e3, ls=":", color="gray", lw=1)
ax.text(uy * 1e3, ax.get_ylim()[1], "  uy", va="top", color="gray")
ax.set_xlabel("u  [mm]"); ax.set_ylabel("F  [kN]")
ax.set_title("Pushover capacity curve — Steel01 perfectly plastic")
ax.legend(); ax.grid(True, alpha=0.3)
fig.tight_layout()

u_fine = np.linspace(0, u_hist[-1], 200) F_closed = np.where(u_fine < uy, E * A * u_fine / L, Fy) fig, ax = plt.subplots(figsize=(6, 4)) ax.plot(u_fine * 1e3, F_closed * 1e-3, "-", color="k", lw=1, label="closed-form bilinear") ax.plot(u_hist * 1e3, F_hist * 1e-3, "o", ms=3, color="C0", label="FEM (capture)") ax.axvline(uy * 1e3, ls=":", color="gray", lw=1) ax.text(uy * 1e3, ax.get_ylim()[1], " uy", va="top", color="gray") ax.set_xlabel("u [mm]"); ax.set_ylabel("F [kN]") ax.set_title("Pushover capacity curve — Steel01 perfectly plastic") ax.legend(); ax.grid(True, alpha=0.3) fig.tight_layout()

9. (Optional) Open the post-solve viewer¶

The viewer's scrubber walks the capacity curve step by step — useful for visualising yield propagation in a richer model than this 1-D bar.

In [ ]:

results.viewer(blocking=False)

results.viewer(blocking=False)

What this unlocks¶

• Multi-step nonlinear analysis through a single spec.capture context. The recorder fires once per converged step; the slab comes back as a (T, N) array ready to plot.
• Capacity curves from slabs. The classical $(u, F)$ pair drops out of two .get() calls — no nodeDisp / nodeReaction accumulation in user code.
• DisplacementControl pattern — a unit reference load defines the controlled DOF; the integrator increments displacement and the analysis returns the reaction.

This closes the curated EOS notebook gallery. The same spec.capture / Results loop scales to richer pushovers — fiber sections, full frames, hysteretic cycles.

In [ ]:

g_ctx.__exit__(None, None, None)

g_ctx.__exit__(None, None, None)