01 — Hello Plate (Plane-Stress Uniaxial Tension)¶

Curriculum slot: Tier 1, slot 01.
Prerequisite: none — this is the entry point.
Strategy used for results: spec.capture(...) — apeGmsh's broadest path.

Problem¶

A square plate of side $L$ and unit thickness is

• fixed in $x$ on the left edge (symmetry),
• anchored in $y$ at the bottom-left corner only (kills rigid-body translation),
• pulled by a uniform horizontal traction $\sigma$ on the right edge.

                 σ [Pa]  →
          ┌─────────────────→┐        +y
          │                  │         │
      ux=0│                  │         └── +x
          │      plate       │
          │      L × L       │
          │   (plane stress) │
          ●──────────────────┘    ← uy = 0 at one corner only

Plane-stress uniaxial tension is a constant-stress field; linear triangles represent it exactly, so the closed-form right-edge displacement is

$$ u_{x,\text{right}} \;=\; \dfrac{\sigma\,L}{E}. $$

Equilibrium gives the total left-edge reaction

$$ R_{x,\text{left}} \;=\; -\,\sigma\,L\,t. $$

What this notebook teaches¶

The complete apeGmsh ↔ openseespy ↔ Results loop on a 2-D continuum:

1. g.model.geometry — build the plate.
2. g.physical — name the boundaries.
3. g.mesh — triangulate.
4. g.mesh.queries.get_fem_data(dim=2) — broker.
5. Vanilla openseespy — nodes / elements / materials / fix / load.
6. apeGmsh recorders — declare what to record by physical-group name.
7. spec.capture(...) — wrap ops.analyze(...), write a native .h5.
8. Results.from_native(...) — read back, verify, plot in-notebook.
9. results.viewer(...) — interactive post-solve view (subprocess).

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 / loading (consistent units: m, Pa, N)
L     = 1.0
thk   = 1.0
sigma = 1.0e6     # right-edge traction

# Material (linear elastic, isotropic)
E  = 2.1e11
nu = 0.3

# Mesh
LC = L / 20.0

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 / loading (consistent units: m, Pa, N) L = 1.0 thk = 1.0 sigma = 1.0e6 # right-edge traction # Material (linear elastic, isotropic) E = 2.1e11 nu = 0.3 # Mesh LC = L / 20.0

2. Geometry¶

Four corner points → four boundary lines → one curve loop → one plane surface. sync() is mandatory before physical groups or meshing can see the geometry.

In [ ]:

g = apeGmsh(model_name="01_hello_plate", verbose=False)
g.begin()

p_BL = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC)
p_BR = g.model.geometry.add_point(L,   0.0, 0.0, lc=LC)
p_TR = g.model.geometry.add_point(L,   L,   0.0, lc=LC)
p_TL = g.model.geometry.add_point(0.0, L,   0.0, lc=LC)

l_B = g.model.geometry.add_line(p_BL, p_BR)    # bottom
l_R = g.model.geometry.add_line(p_BR, p_TR)    # right (loaded)
l_T = g.model.geometry.add_line(p_TR, p_TL)    # top
l_L = g.model.geometry.add_line(p_TL, p_BL)    # left  (fixed in x)

loop    = g.model.geometry.add_curve_loop([l_B, l_R, l_T, l_L])
surface = g.model.geometry.add_plane_surface([loop])

g.model.sync()

g = apeGmsh(model_name="01_hello_plate", verbose=False) g.begin() p_BL = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC) p_BR = g.model.geometry.add_point(L, 0.0, 0.0, lc=LC) p_TR = g.model.geometry.add_point(L, L, 0.0, lc=LC) p_TL = g.model.geometry.add_point(0.0, L, 0.0, lc=LC) l_B = g.model.geometry.add_line(p_BL, p_BR) # bottom l_R = g.model.geometry.add_line(p_BR, p_TR) # right (loaded) l_T = g.model.geometry.add_line(p_TR, p_TL) # top l_L = g.model.geometry.add_line(p_TL, p_BL) # left (fixed in x) loop = g.model.geometry.add_curve_loop([l_B, l_R, l_T, l_L]) surface = g.model.geometry.add_plane_surface([loop]) g.model.sync()

3. Physical groups¶

Naming the boundaries unlocks the rest of the pipeline — the same names we use to apply BCs are the names we'll use to declare recorders and read results back.

In [ ]:

g.physical.add(0, [p_BL],    name="corner")    # uy = 0 anchor
g.physical.add(1, [l_L],     name="left")      # ux = 0 symmetry
g.physical.add(1, [l_R],     name="right")     # σ traction
g.physical.add(2, [surface], name="plate")     # the body

g.physical.add(0, [p_BL], name="corner") # uy = 0 anchor g.physical.add(1, [l_L], name="left") # ux = 0 symmetry g.physical.add(1, [l_R], name="right") # σ traction g.physical.add(2, [surface], name="plate") # the body

4. Mesh¶

In [ ]:

g.mesh.generation.generate(2)
fem = g.mesh.queries.get_fem_data()
print(f"mesh built: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

g.mesh.generation.generate(2) fem = g.mesh.queries.get_fem_data() print(f"mesh built: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

5. Build the OpenSees model — vanilla openseespy¶

The broker hands us clean ID arrays and named groups; from there we drive ops.* directly. No magic — tri31 plane-stress triangles with ElasticIsotropic, fixed left edge, anchored corner, tributary-length right-edge load.

In [ ]:

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

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

# material
MAT_TAG = 1
ops.nDMaterial("ElasticIsotropic", MAT_TAG, E, nu)

# elements (linear plane-stress triangle)
for group in fem.elements.get(target="plate"):
    for eid, conn in zip(group.ids, group.connectivity):
        ops.element(
            "tri31", int(eid),
            int(conn[0]), int(conn[1]), int(conn[2]),
            thk, "PlaneStress", MAT_TAG,
        )

# boundary conditions
for n in fem.nodes.get(target="left").ids:
    ops.fix(int(n), 1, 0)               # ux = 0
for n in fem.nodes.get(target="corner").ids:
    ops.fix(int(n), 0, 1)               # uy = 0 (corner only)

# right-edge traction — distribute by tributary length on y-sorted nodes
right_ids = np.asarray(fem.nodes.get(target="right").ids, dtype=np.int64)
node_y = {int(n): float(fem.nodes.coords[i, 1])
          for i, n in enumerate(fem.nodes.ids)}
y_vals  = np.array([node_y[int(n)] for n in right_ids])
order   = np.argsort(y_vals)
right_sorted = right_ids[order]
y_sorted     = y_vals[order]

seg  = np.diff(y_sorted)
trib = np.empty_like(y_sorted)
trib[0]    = seg[0] / 2.0
trib[-1]   = seg[-1] / 2.0
trib[1:-1] = (seg[:-1] + seg[1:]) / 2.0
assert abs(trib.sum() - L) < 1e-9, f"tributary sum {trib.sum()} != {L}"

ops.timeSeries("Constant", 1)
ops.pattern("Plain", 1, 1)
for nid, t in zip(right_sorted, trib):
    ops.load(int(nid), sigma * float(t) * thk, 0.0)

ops.wipe() ops.model("basic", "-ndm", 2, "-ndf", 2) # nodes for nid, xyz in fem.nodes.get(): ops.node(int(nid), float(xyz[0]), float(xyz[1])) # material MAT_TAG = 1 ops.nDMaterial("ElasticIsotropic", MAT_TAG, E, nu) # elements (linear plane-stress triangle) for group in fem.elements.get(target="plate"): for eid, conn in zip(group.ids, group.connectivity): ops.element( "tri31", int(eid), int(conn[0]), int(conn[1]), int(conn[2]), thk, "PlaneStress", MAT_TAG, ) # boundary conditions for n in fem.nodes.get(target="left").ids: ops.fix(int(n), 1, 0) # ux = 0 for n in fem.nodes.get(target="corner").ids: ops.fix(int(n), 0, 1) # uy = 0 (corner only) # right-edge traction — distribute by tributary length on y-sorted nodes right_ids = np.asarray(fem.nodes.get(target="right").ids, dtype=np.int64) node_y = {int(n): float(fem.nodes.coords[i, 1]) for i, n in enumerate(fem.nodes.ids)} y_vals = np.array([node_y[int(n)] for n in right_ids]) order = np.argsort(y_vals) right_sorted = right_ids[order] y_sorted = y_vals[order] seg = np.diff(y_sorted) trib = np.empty_like(y_sorted) trib[0] = seg[0] / 2.0 trib[-1] = seg[-1] / 2.0 trib[1:-1] = (seg[:-1] + seg[1:]) / 2.0 assert abs(trib.sum() - L) < 1e-9, f"tributary sum {trib.sum()} != {L}" ops.timeSeries("Constant", 1) ops.pattern("Plain", 1, 1) for nid, t in zip(right_sorted, trib): ops.load(int(nid), sigma * float(t) * thk, 0.0)

6. Declare recorders¶

Two named records:

• displacement on all nodes (no pg=) — broad capture so we can plot deformed shapes and post-filter by location.
• reaction_force on pg="left" — targeted, because reactions only matter on the fixed boundary.

Selectors mirror the rest of the library; we instantiate a Recorders() directly and pass ndm / ndf straight to resolve(...) so the resolver knows the DOF space.

In [ ]:

recorders = Recorders()

recorders.nodes(
    components=["displacement"], name="u_all",   # all nodes
)
recorders.nodes(
    components=["reaction_force"], pg="left", name="R_left",
)

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

recorders = Recorders() recorders.nodes( components=["displacement"], name="u_all", # all nodes ) recorders.nodes( components=["reaction_force"], pg="left", name="R_left", ) spec = recorders.resolve(fem, ndm=2, ndf=2) print(spec)

7. Run the analysis with spec.capture(...)¶

Wrap the analysis recipe in a context manager. After each successful ops.analyze(...), cap.step(t=...) queries the live ops domain for everything in the spec. Closing the stage flushes a chunked HDF5 write.

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=2, ndf=2) as cap:
    cap.begin_stage("static", kind="static")

    ops.system("BandGeneral")
    ops.numberer("Plain")
    ops.constraints("Plain")
    ops.test("NormUnbalance", 1e-10, 10)
    ops.algorithm("Linear")
    ops.integrator("LoadControl", 1.0)
    ops.analysis("Static")

    status = ops.analyze(1)
    assert status == 0, f"ops.analyze returned {status}"
    cap.step(t=ops.getTime())

    cap.end_stage()

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=2, ndf=2) as cap: cap.begin_stage("static", kind="static") ops.system("BandGeneral") ops.numberer("Plain") ops.constraints("Plain") ops.test("NormUnbalance", 1e-10, 10) ops.algorithm("Linear") ops.integrator("LoadControl", 1.0) ops.analysis("Static") status = ops.analyze(1) assert status == 0, f"ops.analyze returned {status}" cap.step(t=ops.getTime()) cap.end_stage() print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

8. Read results back — the Results API¶

Open the file, bind it to the same fem, and ask for what we want by physical-group name. No raw HDF5 paths anywhere 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 [ ]:

# Right-edge displacement (verification) — filter post-hoc since
# the recorder captured all nodes.
right_ids = np.asarray(fem.nodes.get(target="right").ids, dtype=np.int64)
u_right = results.nodes.get(
    component="displacement_x",
    ids=right_ids,
    time=-1,
)
u_avg     = float(u_right.values.mean())
u_closed  = sigma * L / E
print(f"u_x at right edge   :  {u_avg:.6e}  m   (computed)")
print(f"σ L / E            :  {u_closed:.6e}  m   (closed form)")
print(f"relative error      :  {abs(u_avg - u_closed) / u_closed:.2e}")

# Right-edge displacement (verification) — filter post-hoc since # the recorder captured all nodes. right_ids = np.asarray(fem.nodes.get(target="right").ids, dtype=np.int64) u_right = results.nodes.get( component="displacement_x", ids=right_ids, time=-1, ) u_avg = float(u_right.values.mean()) u_closed = sigma * L / E print(f"u_x at right edge : {u_avg:.6e} m (computed)") print(f"σ L / E : {u_closed:.6e} m (closed form)") print(f"relative error : {abs(u_avg - u_closed) / u_closed:.2e}")

In [ ]:

# Left-edge reaction (equilibrium)
R_left = results.nodes.get(
    component="reaction_force_x", pg="left", time=-1,
)
Rx_total  = float(R_left.values.sum())
Rx_closed = -sigma * L * thk
print(f"Σ R_x at left edge   :  {Rx_total:.6e}  N   (computed)")
print(f"−σ L t              :  {Rx_closed:.6e}  N   (closed form)")
print(f"relative error      :  {abs(Rx_total - Rx_closed) / abs(Rx_closed):.2e}")

# Left-edge reaction (equilibrium) R_left = results.nodes.get( component="reaction_force_x", pg="left", time=-1, ) Rx_total = float(R_left.values.sum()) Rx_closed = -sigma * L * thk print(f"Σ R_x at left edge : {Rx_total:.6e} N (computed)") print(f"−σ L t : {Rx_closed:.6e} N (closed form)") print(f"relative error : {abs(Rx_total - Rx_closed) / abs(Rx_closed):.2e}")

9. Plot in-notebook¶

Two views: the right-edge displacement profile (uniform along $y$, as uniaxial tension predicts) and the deformed-shape outline.

In [ ]:

# Right-edge u_x along y
y_lookup = dict(zip(
    np.asarray(fem.nodes.ids, dtype=np.int64).tolist(),
    np.asarray(fem.nodes.coords)[:, 1].tolist(),
))
y_at_right = np.array([y_lookup[int(n)] for n in u_right.node_ids])
u_at_right = u_right.values[0]
order      = np.argsort(y_at_right)

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(u_at_right[order], y_at_right[order], "o-", ms=4, label="computed")
ax.axvline(u_closed, ls="--", color="k", lw=1, label="closed form  σL/E")
ax.set_xlabel("u_x  [m]"); ax.set_ylabel("y  [m]")
ax.set_title("Right-edge displacement (uniform along y)")
ax.legend(); fig.tight_layout()

# Right-edge u_x along y y_lookup = dict(zip( np.asarray(fem.nodes.ids, dtype=np.int64).tolist(), np.asarray(fem.nodes.coords)[:, 1].tolist(), )) y_at_right = np.array([y_lookup[int(n)] for n in u_right.node_ids]) u_at_right = u_right.values[0] order = np.argsort(y_at_right) fig, ax = plt.subplots(figsize=(5, 4)) ax.plot(u_at_right[order], y_at_right[order], "o-", ms=4, label="computed") ax.axvline(u_closed, ls="--", color="k", lw=1, label="closed form σL/E") ax.set_xlabel("u_x [m]"); ax.set_ylabel("y [m]") ax.set_title("Right-edge displacement (uniform along y)") ax.legend(); fig.tight_layout()

In [ ]:

# Reference vs deformed (boundary scatter, scaled)
scale = 5.0e3
ux_all = results.nodes.get(component="displacement_x", time=-1).values[0]
uy_all = results.nodes.get(component="displacement_y", time=-1).values[0]
X = fem.nodes.coords[:, 0] + scale * ux_all
Y = fem.nodes.coords[:, 1] + scale * uy_all

fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(fem.nodes.coords[:, 0], fem.nodes.coords[:, 1],
           s=2, color="lightgray", label="reference")
ax.scatter(X, Y, s=4, color="C1", label=f"deformed (×{int(scale):d})")
ax.set_aspect("equal")
ax.set_xlabel("x  [m]"); ax.set_ylabel("y  [m]")
ax.set_title("Reference vs deformed (uniaxial tension)")
ax.legend(); fig.tight_layout()

# Reference vs deformed (boundary scatter, scaled) scale = 5.0e3 ux_all = results.nodes.get(component="displacement_x", time=-1).values[0] uy_all = results.nodes.get(component="displacement_y", time=-1).values[0] X = fem.nodes.coords[:, 0] + scale * ux_all Y = fem.nodes.coords[:, 1] + scale * uy_all fig, ax = plt.subplots(figsize=(5, 5)) ax.scatter(fem.nodes.coords[:, 0], fem.nodes.coords[:, 1], s=2, color="lightgray", label="reference") ax.scatter(X, Y, s=4, color="C1", label=f"deformed (×{int(scale):d})") ax.set_aspect("equal") ax.set_xlabel("x [m]"); ax.set_ylabel("y [m]") ax.set_title("Reference vs deformed (uniaxial tension)") ax.legend(); fig.tight_layout()

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

Subprocess mode — the notebook stays alive while the viewer window is open. Set APEGMSH_SKIP_VIEWER=1 to skip in headless / CI runs.

In [ ]:

results.viewer(blocking=False)

results.viewer(blocking=False)

What this unlocks¶

You've walked the full pipeline:

• Build geometry / mesh / FEM with apeGmsh.
• Drive OpenSees with vanilla ops.* calls.
• Declare recorders by physical-group name.
• Capture the analysis with spec.capture(...).
• Read results via Results.nodes.get(...) — same selector vocabulary as the declaration side.
• Plot / view in-notebook from the slabs without touching ops.nodeDisp in user code.

The same template scales to bigger meshes, more recorders, multi-stage (gravity → dynamic), modal analysis (cap.capture_modes(...)), and the read-side spatial filters (results.nodes.in_box(...)).

In [ ]:

g.end()

g.end()