04 — 2D Portal Frame (Lateral Load)¶

Curriculum slot: Tier 1, slot 04.
Prerequisite: 02 — Cantilever Beam.
Strategy used for results: spec.capture(...) — apeGmsh's broadest path.

Problem¶

A single-bay planar portal frame: two columns of height $h$ fixed at the base, connected at the top by a horizontal beam of length $b$. The beam is modelled effectively rigid (its stiffness is bumped well above the columns) so the classical "rigid-beam portal" closed form applies. A lateral load $H$ is applied at the top-left joint.

   H→  ┌────────────────────────┐      +z
       │                        │       │
       │                        │       └── +x
    h  │      (columns)         │
       │                        │
       │                        │
       ●                        ●
   ────────  fixed bases  ─────────
              b = 4 m

Each column behaves as a fixed-fixed shear element with lateral stiffness $k_{\text{col}} = 12 E I_c / h^3$. The two columns act in parallel ($K = 2 k_{\text{col}}$), so the classical drift is

$$ \Delta \;=\; \dfrac{H\,h^3}{24\,E\,I_c}, $$

and each column carries a base moment of $H h / 4$.

Why this notebook is interesting¶

• Multi-element model with two assemblies of different cross-section — columns and beam — driven from named physical groups.
• First example where the FEM and the closed-form check don't match to round-off. The classical formula assumes axial-rigid columns and a point joint; elasticBeamColumn honours neither. The FEM is more correct, and the discrepancy is the lesson (see Section 8.5).

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 / load (consistent units: m, Pa, N)
h = 3.0          # column height
b = 4.0          # beam span
H = 10_000.0     # lateral load at top-left joint

# Material
E  = 2.1e11
nu = 0.3
G  = E / (2.0 * (1.0 + nu))

# Column cross-section
A_col  = 1.0e-3
Iy_col = 1.0e-5
Iz_col = 1.0e-5     # in-plane bending — controls the frame drift
J_col  = 2.0e-5

# Beam cross-section — bumped large so it acts effectively rigid.
A_beam, Iy_beam, Iz_beam, J_beam = 1.0, 1.0, 1.0, 2.0

# Mesh density
LC = min(h, b) / 10.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 / load (consistent units: m, Pa, N) h = 3.0 # column height b = 4.0 # beam span H = 10_000.0 # lateral load at top-left joint # Material E = 2.1e11 nu = 0.3 G = E / (2.0 * (1.0 + nu)) # Column cross-section A_col = 1.0e-3 Iy_col = 1.0e-5 Iz_col = 1.0e-5 # in-plane bending — controls the frame drift J_col = 2.0e-5 # Beam cross-section — bumped large so it acts effectively rigid. A_beam, Iy_beam, Iz_beam, J_beam = 1.0, 1.0, 1.0, 2.0 # Mesh density LC = min(h, b) / 10.0

2. Geometry¶

Four corner points → two column lines + one beam line. The frame lives in the $xz$-plane; the second axis $y$ is unused / restrained later.

In [ ]:

g = apeGmsh(model_name="04_portal_frame_2D", verbose=False)
g.begin()

p_BL = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC)   # base left
p_BR = g.model.geometry.add_point(b,   0.0, 0.0, lc=LC)   # base right
p_TL = g.model.geometry.add_point(0.0, 0.0, h,   lc=LC)   # top-left (loaded)
p_TR = g.model.geometry.add_point(b,   0.0, h,   lc=LC)   # top-right

col_L = g.model.geometry.add_line(p_BL, p_TL)
col_R = g.model.geometry.add_line(p_BR, p_TR)
beam  = g.model.geometry.add_line(p_TL, p_TR)

g.model.sync()

g = apeGmsh(model_name="04_portal_frame_2D", verbose=False) g.begin() p_BL = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC) # base left p_BR = g.model.geometry.add_point(b, 0.0, 0.0, lc=LC) # base right p_TL = g.model.geometry.add_point(0.0, 0.0, h, lc=LC) # top-left (loaded) p_TR = g.model.geometry.add_point(b, 0.0, h, lc=LC) # top-right col_L = g.model.geometry.add_line(p_BL, p_TL) col_R = g.model.geometry.add_line(p_BR, p_TR) beam = g.model.geometry.add_line(p_TL, p_TR) g.model.sync()

3. Physical groups¶

Five names: bases (two fixed supports), top_left (loaded), top_right, columns (line group for the two columns), beam.

In [ ]:

g.physical.add(0, [p_BL, p_BR],   name="bases")
g.physical.add(0, [p_TL],         name="top_left")
g.physical.add(0, [p_TR],         name="top_right")
g.physical.add(1, [col_L, col_R], name="columns")
g.physical.add(1, [beam],         name="beam")

g.physical.add(0, [p_BL, p_BR], name="bases") g.physical.add(0, [p_TL], name="top_left") g.physical.add(0, [p_TR], name="top_right") g.physical.add(1, [col_L, col_R], name="columns") g.physical.add(1, [beam], name="beam")

4. Mesh¶

In [ ]:

g.mesh.generation.generate(1)
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(1) 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¶

Two element groups (columns + beam) get different sections via separate fem.elements.get(target=...) loops. The frame lives in the $xz$-plane, so non-base nodes get the out-of-plane DOFs (uy, rx, rz) restrained — otherwise the frame has zero-energy out-of-plane rigid-body modes.

In [ ]:

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

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

# geometric transformations — same orientation works for both groups here
TRANSF_COL  = 1
TRANSF_BEAM = 2
ops.geomTransf("Linear", TRANSF_COL,  0.0, 1.0, 0.0)
ops.geomTransf("Linear", TRANSF_BEAM, 0.0, 1.0, 0.0)

# columns (smaller section)
for group in fem.elements.get(target="columns"):
    for eid, conn in zip(group.ids, group.connectivity):
        ops.element(
            "elasticBeamColumn", int(eid),
            int(conn[0]), int(conn[1]),
            A_col, E, G, J_col, Iy_col, Iz_col, TRANSF_COL,
        )

# beam (effectively rigid section)
for group in fem.elements.get(target="beam"):
    for eid, conn in zip(group.ids, group.connectivity):
        ops.element(
            "elasticBeamColumn", int(eid),
            int(conn[0]), int(conn[1]),
            A_beam, E, G, J_beam, Iy_beam, Iz_beam, TRANSF_BEAM,
        )

# fully fix the bases
base_set = {int(n) for n in fem.nodes.get(target="bases").ids}
for n in base_set:
    ops.fix(int(n), 1, 1, 1, 1, 1, 1)

# planar restraint on every other node (uy, rx, rz fixed)
for nid, _ in fem.nodes.get():
    if int(nid) in base_set:
        continue
    ops.fix(int(nid), 0, 1, 0, 1, 0, 1)

# lateral point load at the top-left joint (+x)
ops.timeSeries("Constant", 1)
ops.pattern("Plain", 1, 1)
for n in fem.nodes.get(target="top_left").ids:
    ops.load(int(n), H, 0.0, 0.0, 0.0, 0.0, 0.0)

ops.wipe() ops.model("basic", "-ndm", 3, "-ndf", 6) # nodes for nid, xyz in fem.nodes.get(): ops.node(int(nid), float(xyz[0]), float(xyz[1]), float(xyz[2])) # geometric transformations — same orientation works for both groups here TRANSF_COL = 1 TRANSF_BEAM = 2 ops.geomTransf("Linear", TRANSF_COL, 0.0, 1.0, 0.0) ops.geomTransf("Linear", TRANSF_BEAM, 0.0, 1.0, 0.0) # columns (smaller section) for group in fem.elements.get(target="columns"): for eid, conn in zip(group.ids, group.connectivity): ops.element( "elasticBeamColumn", int(eid), int(conn[0]), int(conn[1]), A_col, E, G, J_col, Iy_col, Iz_col, TRANSF_COL, ) # beam (effectively rigid section) for group in fem.elements.get(target="beam"): for eid, conn in zip(group.ids, group.connectivity): ops.element( "elasticBeamColumn", int(eid), int(conn[0]), int(conn[1]), A_beam, E, G, J_beam, Iy_beam, Iz_beam, TRANSF_BEAM, ) # fully fix the bases base_set = {int(n) for n in fem.nodes.get(target="bases").ids} for n in base_set: ops.fix(int(n), 1, 1, 1, 1, 1, 1) # planar restraint on every other node (uy, rx, rz fixed) for nid, _ in fem.nodes.get(): if int(nid) in base_set: continue ops.fix(int(nid), 0, 1, 0, 1, 0, 1) # lateral point load at the top-left joint (+x) ops.timeSeries("Constant", 1) ops.pattern("Plain", 1, 1) for n in fem.nodes.get(target="top_left").ids: ops.load(int(n), H, 0.0, 0.0, 0.0, 0.0, 0.0)

6. Declare recorders¶

• displacement on all nodes — for the deformed-shape plot.
• reaction on pg="bases" — base shears + moments per column.

In [ ]:

recorders = Recorders()

recorders.nodes(
    components=["displacement"], name="u_all",
)
recorders.nodes(
    components=["reaction"], pg="bases", name="R_bases",
)

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

recorders = Recorders() recorders.nodes( components=["displacement"], name="u_all", ) recorders.nodes( components=["reaction"], pg="bases", name="R_bases", ) spec = recorders.resolve(fem, ndm=3, ndf=6) print(spec)

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

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=3, ndf=6) 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=3, ndf=6) 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 — verifications¶

Two things to check:

1. Lateral drift at the loaded top-left joint vs $\Delta = Hh^3/(24EI_c)$.
2. Base moment at one column base vs $M_{\text{base}} = Hh/4$.

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 [ ]:

# Lateral drift at top_left
tl_id = int(next(iter(fem.nodes.get(target="top_left").ids)))
drift_fem    = float(results.nodes.get(
    component="displacement_x", ids=[tl_id], time=-1,
).values[0, 0])
drift_closed = H * h**3 / (24.0 * E * Iz_col)
err_drift    = abs(drift_fem - drift_closed) / drift_closed * 100.0

print("Top-left lateral drift")
print(f"  FEM                 :  {drift_fem:.6e}  m")
print(f"  H h^3 / (24 E I_c)  :  {drift_closed:.6e}  m   (closed form)")
print(f"  Error               :  {err_drift:.4f} %")

# Lateral drift at top_left tl_id = int(next(iter(fem.nodes.get(target="top_left").ids))) drift_fem = float(results.nodes.get( component="displacement_x", ids=[tl_id], time=-1, ).values[0, 0]) drift_closed = H * h**3 / (24.0 * E * Iz_col) err_drift = abs(drift_fem - drift_closed) / drift_closed * 100.0 print("Top-left lateral drift") print(f" FEM : {drift_fem:.6e} m") print(f" H h^3 / (24 E I_c) : {drift_closed:.6e} m (closed form)") print(f" Error : {err_drift:.4f} %")

In [ ]:

# Base moment per column — magnitude of M_y reaction at each base
M_y = results.nodes.get(
    component="reaction_moment_y", pg="bases", time=-1,
)
M_per_col = float(np.abs(M_y.values[0]).max())
M_closed  = H * h / 4.0
err_M     = abs(M_per_col - M_closed) / M_closed * 100.0

print("Column base moment (magnitude per column)")
print(f"  FEM         :  {M_per_col:.6e}  N·m")
print(f"  H h / 4     :  {M_closed:.6e}  N·m   (closed form)")
print(f"  Error       :  {err_M:.4f} %")

# Base moment per column — magnitude of M_y reaction at each base M_y = results.nodes.get( component="reaction_moment_y", pg="bases", time=-1, ) M_per_col = float(np.abs(M_y.values[0]).max()) M_closed = H * h / 4.0 err_M = abs(M_per_col - M_closed) / M_closed * 100.0 print("Column base moment (magnitude per column)") print(f" FEM : {M_per_col:.6e} N·m") print(f" H h / 4 : {M_closed:.6e} N·m (closed form)") print(f" Error : {err_M:.4f} %")

8.5 Why ~0.75% drift error instead of 0?¶

Slot 02 matched analytical theory to machine precision. This one doesn't — the residual ~0.75% is inherent, not a meshing or tolerance issue. Bumping the beam stiffness even higher leaves the error unchanged.

The classical formula $\Delta = H h^3 / (24 E I_c)$ assumes:

1. Rigid beam — top nodes translate together with no relative rotation (satisfied by our 10⁵× beam).
2. Columns have flexural DOFs only — no axial deformation, no shear.
3. Connections are points — no joint-zone effect.

elasticBeamColumn honours none of these — the columns compress and extend under the overturning axial couple, the Timoshenko-style stiffness adds shear flex, and the joint participates as a regular 6-DOF node. The FEM drift is therefore larger than the classical formula, and closer to reality. This is the first curriculum example where the FEM is more correct than the analytical check.

9. Plot the deformed frame¶

Pull displacement_x and displacement_z at every node and overlay the deformed positions on the reference geometry.

In [ ]:

ux_all = results.nodes.get(component="displacement_x", time=-1).values[0]
uz_all = results.nodes.get(component="displacement_z", time=-1).values[0]

scale = 30.0          # visualization scale
X_ref = fem.nodes.coords[:, 0]
Z_ref = fem.nodes.coords[:, 2]
X_def = X_ref + scale * ux_all
Z_def = Z_ref + scale * uz_all

fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(X_ref, Z_ref, s=8, color="lightgray", label="reference")
ax.scatter(X_def, Z_def, s=12, color="C1", label=f"deformed (×{scale:g})")
ax.axhline(0.0, ls=":", color="gray", lw=1)
ax.set_aspect("equal")
ax.set_xlabel("x  [m]"); ax.set_ylabel("z  [m]")
ax.set_title("Portal frame — reference vs deformed (lateral load)")
ax.legend()
fig.tight_layout()

ux_all = results.nodes.get(component="displacement_x", time=-1).values[0] uz_all = results.nodes.get(component="displacement_z", time=-1).values[0] scale = 30.0 # visualization scale X_ref = fem.nodes.coords[:, 0] Z_ref = fem.nodes.coords[:, 2] X_def = X_ref + scale * ux_all Z_def = Z_ref + scale * uz_all fig, ax = plt.subplots(figsize=(6, 5)) ax.scatter(X_ref, Z_ref, s=8, color="lightgray", label="reference") ax.scatter(X_def, Z_def, s=12, color="C1", label=f"deformed (×{scale:g})") ax.axhline(0.0, ls=":", color="gray", lw=1) ax.set_aspect("equal") ax.set_xlabel("x [m]"); ax.set_ylabel("z [m]") ax.set_title("Portal frame — reference vs deformed (lateral load)") ax.legend() fig.tight_layout()

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

In [ ]:

results.viewer(blocking=False)

results.viewer(blocking=False)

What this unlocks¶

• Multi-element model with multiple sections assigned through named physical groups — columns get one cross-section, beam another. Same idiom scales to dozens of element groups.
• Two independent geomTransf slots — the minimum needed for any non-collinear beam network.
• Out-of-plane restraint idiom for 2D-in-3D analyses (fix(n, 0, 1, 0, 1, 0, 1) on every non-base node).
• First example where the FEM is closer to reality than a textbook closed form. Useful pattern for sanity-checking nonlinear models against linear-elastic baselines downstream.

Next: notebook 05 — labels and physical groups, the naming patterns the rest of the curriculum builds on.

In [ ]:

g.end()

g.end()