12 — Interface Coupling via equal_dof¶

Curriculum slot: Tier 3, slot 12.
Prerequisite: 10 — Parts Basics, 11 — Boolean Assembly.
Strategy used for results: spec.capture(...).

Purpose¶

Slot 11 showed how a boolean fragment produces a conformal mesh at the interface of two parts — the CAD itself becomes one body, so Gmsh places a single shared node along the common edge.

This slot covers the alternative: the two parts stay geometrically separate. Each part meshes independently, the junction has two coincident-but-distinct mesh nodes, and a constraint glues their DOFs together.

The simplest form is equal_dof — coordinate-matched nodes share all requested DOFs. For non-matching meshes with no coincident nodes, a tie with shape-function interpolation is needed instead (slot 16).

Problem¶

Two linear-elastic cantilever beams meeting at $x = L_A$. Beam A spans $[0, L_A]$ with its fixed base at $x = 0$. Beam B spans $[L_A, L_A + L_B]$ with a downward tip load $P$. Each beam is a separate part, so the junction is two coincident points ($p_{A1}$ and $p_{B0}$) and after meshing two coincident nodes.

With the junction tied (equal_dof on all 6 DOFs) the structure behaves as a single cantilever of total length $L = L_A + L_B$:

$$ \delta_{\text{tip}} \;=\; -\,\dfrac{P\,L^3}{3\,E\,I}, \qquad L = L_A + L_B. $$

Why this notebook is interesting¶

• g.constraints.equal_dof resolves into per-pair constraint records on fem.nodes.constraints.equal_dofs() — clean read out of the broker and clean emission via ops.equalDOF(...).
• spec.capture records both bases' reactions and the full displacement field — the read side then verifies that the master and slave nodes at the junction have identical displacements (the equal_dof contract).
• Per-part deflection plot via label="beamA" / label="beamB" on the read side.

1. Imports and parameters¶

In [ ]:

from pathlib import Path

import gmsh
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

L_A, L_B = 1.5, 1.5
L = L_A + L_B
P = 10_000.0
E  = 2.1e11
nu = 0.3
G  = E / (2.0 * (1.0 + nu))
A, Iy, Iz, J = 1.0e-3, 1.0e-5, 1.0e-5, 2.0e-5
LC = L / 20.0

from pathlib import Path import gmsh 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 L_A, L_B = 1.5, 1.5 L = L_A + L_B P = 10_000.0 E = 2.1e11 nu = 0.3 G = E / (2.0 * (1.0 + nu)) A, Iy, Iz, J = 1.0e-3, 1.0e-5, 1.0e-5, 2.0e-5 LC = L / 20.0

2. Geometry — two independent beams¶

Each beam is a separate Part — geometrically coincident at $x = L_A$ but topologically distinct. No fragment, no fuse.

In [ ]:

g_ctx = apeGmsh(model_name="12_interface_tie", verbose=False)
g = g_ctx.__enter__()

# Beam A
pA0 = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC)
pA1 = g.model.geometry.add_point(L_A, 0.0, 0.0, lc=LC)     # junction, side A
lnA = g.model.geometry.add_line(pA0, pA1)

# Beam B — a DISTINCT point at the junction
pB0 = g.model.geometry.add_point(L_A, 0.0, 0.0, lc=LC)     # junction, side B
pB1 = g.model.geometry.add_point(L,   0.0, 0.0, lc=LC)
lnB = g.model.geometry.add_line(pB0, pB1)

g.model.sync()

# Register the two parts (no fragment — keep bodies distinct)
instA = g.parts.register("beamA", [(1, lnA), (0, pA0), (0, pA1)])
instB = g.parts.register("beamB", [(1, lnB), (0, pB0), (0, pB1)])

g_ctx = apeGmsh(model_name="12_interface_tie", verbose=False) g = g_ctx.__enter__() # Beam A pA0 = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC) pA1 = g.model.geometry.add_point(L_A, 0.0, 0.0, lc=LC) # junction, side A lnA = g.model.geometry.add_line(pA0, pA1) # Beam B — a DISTINCT point at the junction pB0 = g.model.geometry.add_point(L_A, 0.0, 0.0, lc=LC) # junction, side B pB1 = g.model.geometry.add_point(L, 0.0, 0.0, lc=LC) lnB = g.model.geometry.add_line(pB0, pB1) g.model.sync() # Register the two parts (no fragment — keep bodies distinct) instA = g.parts.register("beamA", [(1, lnA), (0, pA0), (0, pA1)]) instB = g.parts.register("beamB", [(1, lnB), (0, pB0), (0, pB1)])

3. Physical groups + the tie¶

In [ ]:

g.physical.add(0, [pA0], name="base")
g.physical.add(0, [pB1], name="tip")

g.physical.add(0, [pA0], name="base") g.physical.add(0, [pB1], name="tip")

4. Declare the interface tie¶

g.constraints.equal_dof matches coincident nodes across the two parts' interface entities and emits one NodePairRecord per match. Bounding the search to specific entities prevents accidental matches elsewhere in the model.

In [ ]:

g.constraints.equal_dof(
    master_label="beamA",
    slave_label="beamB",
    master_entities=[(0, pA1)],       # junction point on side A
    slave_entities=[(0, pB0)],        # junction point on side B
    tolerance=1e-9,
    dofs=[1, 2, 3, 4, 5, 6],          # couple all 6 DOFs
)

g.constraints.equal_dof( master_label="beamA", slave_label="beamB", master_entities=[(0, pA1)], # junction point on side A slave_entities=[(0, pB0)], # junction point on side B tolerance=1e-9, dofs=[1, 2, 3, 4, 5, 6], # couple all 6 DOFs )

5. Mesh + inspect resolved constraint records¶

In [ ]:

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")

equal_dof_pairs = list(fem.nodes.constraints.equal_dofs())
print(f"equal_dof pairs resolved: {len(equal_dof_pairs)}")
for rec in equal_dof_pairs:
    print(f"  master {rec.master_node} ↔ slave {rec.slave_node}  dofs={rec.dofs}")

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") equal_dof_pairs = list(fem.nodes.constraints.equal_dofs()) print(f"equal_dof pairs resolved: {len(equal_dof_pairs)}") for rec in equal_dof_pairs: print(f" master {rec.master_node} ↔ slave {rec.slave_node} dofs={rec.dofs}")

6. Build the OpenSees model — vanilla openseespy¶

Standard ingest plus one new step: iterate fem.nodes.constraints.equal_dofs() and emit ops.equalDOF per pair — no manual node-ID matching, no manual coincidence search.

In [ ]:

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

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

ops.geomTransf("Linear", 1, 0.0, 1.0, 0.0)

def line_elements_of(inst):
    out = []
    for tag in inst.entities.get(1, []):
        etypes, etags, enodes = gmsh.model.mesh.getElements(1, int(tag))
        for etype, elist, nlist in zip(etypes, etags, enodes):
            if int(etype) != 1:
                continue
            arr = np.asarray(nlist, dtype=np.int64).reshape(-1, 2)
            for eid, row in zip(elist, arr):
                out.append((int(eid), (int(row[0]), int(row[1]))))
    return out

for inst in (instA, instB):
    for eid, (ni, nj) in line_elements_of(inst):
        ops.element(
            "elasticBeamColumn", eid, ni, nj,
            A, E, G, J, Iy, Iz, 1,
        )

# The interface tie — one ops.equalDOF per resolved pair
for rec in fem.nodes.constraints.equal_dofs():
    ops.equalDOF(int(rec.master_node), int(rec.slave_node), *rec.dofs)

# BCs + tip load
for n in fem.nodes.get(target="base").ids:
    ops.fix(int(n), 1, 1, 1, 1, 1, 1)

ops.timeSeries("Constant", 1)
ops.pattern("Plain", 1, 1)
for n in fem.nodes.get(target="tip").ids:
    ops.load(int(n), 0.0, 0.0, -P, 0.0, 0.0, 0.0)

ops.wipe() ops.model("basic", "-ndm", 3, "-ndf", 6) for nid, xyz in fem.nodes.get(): ops.node(int(nid), float(xyz[0]), float(xyz[1]), float(xyz[2])) ops.geomTransf("Linear", 1, 0.0, 1.0, 0.0) def line_elements_of(inst): out = [] for tag in inst.entities.get(1, []): etypes, etags, enodes = gmsh.model.mesh.getElements(1, int(tag)) for etype, elist, nlist in zip(etypes, etags, enodes): if int(etype) != 1: continue arr = np.asarray(nlist, dtype=np.int64).reshape(-1, 2) for eid, row in zip(elist, arr): out.append((int(eid), (int(row[0]), int(row[1])))) return out for inst in (instA, instB): for eid, (ni, nj) in line_elements_of(inst): ops.element( "elasticBeamColumn", eid, ni, nj, A, E, G, J, Iy, Iz, 1, ) # The interface tie — one ops.equalDOF per resolved pair for rec in fem.nodes.constraints.equal_dofs(): ops.equalDOF(int(rec.master_node), int(rec.slave_node), *rec.dofs) # BCs + tip load for n in fem.nodes.get(target="base").ids: ops.fix(int(n), 1, 1, 1, 1, 1, 1) ops.timeSeries("Constant", 1) ops.pattern("Plain", 1, 1) for n in fem.nodes.get(target="tip").ids: ops.load(int(n), 0.0, 0.0, -P, 0.0, 0.0, 0.0)

7. Declare recorders¶

• displacement on all nodes — for the deflection plot and to cross-check the master/slave continuity at the junction.
• reaction at the base PG — equilibrium check.

In [ ]:

recorders = Recorders()

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

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

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

8. 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")
    ops.analyze(1)
    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") ops.analyze(1) cap.step(t=ops.getTime()) cap.end_stage() print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

9. Read results back — three checks¶

1. Tip deflection vs the single-cantilever closed form.
2. Master/slave displacement equality at the junction — the equal_dof contract (this is the reason for the constraint).
3. Base reaction — equilibrium.

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

# 1. Tip deflection
tip_id = int(next(iter(fem.nodes.get(target="tip").ids)))
u_tip = float(results.nodes.get(
    component="displacement_z", ids=[tip_id], time=-1,
).values[0, 0])
u_closed = -P * L**3 / (3.0 * E * Iz)
print(f"u_z at tip                  :  {u_tip:+.6e}  m")
print(f"−P L^3 / (3 E I)            :  {u_closed:+.6e}  m   (closed form, L=L_A+L_B)")
print(f"relative error              :  {abs(u_tip - u_closed)/abs(u_closed):.2e}")

# 1. Tip deflection tip_id = int(next(iter(fem.nodes.get(target="tip").ids))) u_tip = float(results.nodes.get( component="displacement_z", ids=[tip_id], time=-1, ).values[0, 0]) u_closed = -P * L**3 / (3.0 * E * Iz) print(f"u_z at tip : {u_tip:+.6e} m") print(f"−P L^3 / (3 E I) : {u_closed:+.6e} m (closed form, L=L_A+L_B)") print(f"relative error : {abs(u_tip - u_closed)/abs(u_closed):.2e}")

In [ ]:

# 2. Master / slave equality at the junction
master_id = int(equal_dof_pairs[0].master_node)
slave_id  = int(equal_dof_pairs[0].slave_node)
u_master = float(results.nodes.get(
    component="displacement_z", ids=[master_id], time=-1,
).values[0, 0])
u_slave = float(results.nodes.get(
    component="displacement_z", ids=[slave_id], time=-1,
).values[0, 0])
print(f"u_z at master ({master_id})  :  {u_master:+.6e}  m")
print(f"u_z at slave  ({slave_id})   :  {u_slave:+.6e}  m")
print(f"|u_master − u_slave|        :  {abs(u_master - u_slave):.2e}  m   (should be ~ 0)")

# 2. Master / slave equality at the junction master_id = int(equal_dof_pairs[0].master_node) slave_id = int(equal_dof_pairs[0].slave_node) u_master = float(results.nodes.get( component="displacement_z", ids=[master_id], time=-1, ).values[0, 0]) u_slave = float(results.nodes.get( component="displacement_z", ids=[slave_id], time=-1, ).values[0, 0]) print(f"u_z at master ({master_id}) : {u_master:+.6e} m") print(f"u_z at slave ({slave_id}) : {u_slave:+.6e} m") print(f"|u_master − u_slave| : {abs(u_master - u_slave):.2e} m (should be ~ 0)")

In [ ]:

# 3. Base reaction — equilibrium
R_z = float(results.nodes.get(
    component="reaction_force_z", pg="base", time=-1,
).values.sum())
print(f"Σ R_z at base               :  {R_z:+.6e}  N    (equilibrium: +P = {P:+.6e})")

# 3. Base reaction — equilibrium R_z = float(results.nodes.get( component="reaction_force_z", pg="base", time=-1, ).values.sum()) print(f"Σ R_z at base : {R_z:+.6e} N (equilibrium: +P = {P:+.6e})")

10. Plot the deflected shape — partitioned per part¶

Pull displacement_z for each part separately via label= so the two beams can be drawn in different colours, with the junction highlighted.

In [ ]:

tag_to_idx = {int(t): i for i, t in enumerate(fem.nodes.ids)}

fig, ax = plt.subplots(figsize=(7, 4))
x_fine = np.linspace(0, L, 200)
w_closed = -P * x_fine**2 / (6.0 * E * Iz) * (3.0 * L - x_fine)
ax.plot(x_fine, w_closed, "-", color="k", lw=1, label="closed form")

for i, lbl in enumerate(("beamA", "beamB")):
    # Part labels are FEM-side selectors (resolve through fem.nodes.get(target=lbl)).
    # The read side accepts pg/label/selection/ids — so look up IDs first.
    part_ids = np.asarray(fem.nodes.get(target=lbl).ids, dtype=np.int64)
    slab = results.nodes.get(
        component="displacement_z", ids=part_ids, time=-1,
    )
    x = np.array([fem.nodes.coords[tag_to_idx[int(n)], 0]
                  for n in slab.node_ids])
    u = slab.values[0]
    order = np.argsort(x)
    ax.plot(x[order], u[order], "o-", ms=4, color=f"C{i}", label=lbl)

# Highlight the junction
ax.axvline(L_A, color="gray", ls=":", lw=1)
ax.text(L_A, ax.get_ylim()[1], "  junction (equal_dof tie)", va="top", color="gray")

ax.set_xlabel("x  [m]"); ax.set_ylabel("u_z(x)  [m]")
ax.set_title("Two beams tied at x = L_A — same answer as one cantilever")
ax.legend(); ax.grid(True, alpha=0.3)
fig.tight_layout()

tag_to_idx = {int(t): i for i, t in enumerate(fem.nodes.ids)} fig, ax = plt.subplots(figsize=(7, 4)) x_fine = np.linspace(0, L, 200) w_closed = -P * x_fine**2 / (6.0 * E * Iz) * (3.0 * L - x_fine) ax.plot(x_fine, w_closed, "-", color="k", lw=1, label="closed form") for i, lbl in enumerate(("beamA", "beamB")): # Part labels are FEM-side selectors (resolve through fem.nodes.get(target=lbl)). # The read side accepts pg/label/selection/ids — so look up IDs first. part_ids = np.asarray(fem.nodes.get(target=lbl).ids, dtype=np.int64) slab = results.nodes.get( component="displacement_z", ids=part_ids, time=-1, ) x = np.array([fem.nodes.coords[tag_to_idx[int(n)], 0] for n in slab.node_ids]) u = slab.values[0] order = np.argsort(x) ax.plot(x[order], u[order], "o-", ms=4, color=f"C{i}", label=lbl) # Highlight the junction ax.axvline(L_A, color="gray", ls=":", lw=1) ax.text(L_A, ax.get_ylim()[1], " junction (equal_dof tie)", va="top", color="gray") ax.set_xlabel("x [m]"); ax.set_ylabel("u_z(x) [m]") ax.set_title("Two beams tied at x = L_A — same answer as one cantilever") ax.legend(); ax.grid(True, alpha=0.3) fig.tight_layout()

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

In [ ]:

results.viewer(blocking=False)

results.viewer(blocking=False)

What this unlocks¶

• Interface coupling mechanism — g.constraints.equal_dof finds coincident nodes across parts and produces clean NodePairRecord entries. Emission to OpenSees is one line per pair.
• Unfragmented assembly workflow — keep parts geometrically separate, mesh each independently, tie them through constraints. Preferred whenever the two parts use different element types or material laws.
• Read-side verification of the constraint contract — pulling master and slave displacements through the unified Results API shows directly that equal_dof is doing what it says.
• Scaffolding for slot 16. That slot replaces equal_dof with g.constraints.tie for non-matching meshes, where each slave node projects onto a master element face with shape-function weights. Same resolve / record / ingest chain — only the emission swaps ops.equalDOF for ASDEmbeddedNodeElement.

Next: notebook 17 — modal analysis (capture-only territory).

In [ ]:

g_ctx.__exit__(None, None, None)

g_ctx.__exit__(None, None, None)