Constraints

This page explains how apeGmsh couples regions of a model that don't share nodes — which constraint fits which coupling situation, when declarations turn into node-level records, and why you never hand-write a single equalDOF for the solver.

One equation, four situations

Every multi-point constraint in the library is ultimately the same statement: some slave DOFs are a linear function of some master DOFs, u_slave = C · u_master. What distinguishes the dozen factory methods on g.constraints is how the coefficient matrix gets built, and that follows directly from the geometry of the interface:

• the two sides share node positions → pair co-located nodes directly, or drive a node group rigidly from one master;
• the two sides meet at a surface but the meshes don't match → project slave nodes onto master faces and interpolate through shape functions;
• one part is buried inside another (rebar in concrete) → locate each embedded node inside a host element and couple it to the host's corners;
• the two sides may touch, slide, or separate → a genuine contact formulation, not a kinematic constraint at all.

Pick the situation and the method follows. Everything below walks the four in that order, but first the timing — because it's the same for all of them, and it's the part that makes constraints feel different from other FEM tools.

Declare, then resolve, then emit

Constraints follow the declare-then-resolve rhythm from the mental model. You declare a constraint before the mesh exists, against names — part labels or physical groups, never node tags:

# Stage 1 — declare (pre-mesh)
g.constraints.equal_dof("slab", "column_top", dofs=[1, 2, 3])

# Stage 2 — resolve (happens inside get_fem_data)
fem = g.mesh.queries.get_fem_data(dim=3)

# Resolved records are now on the broker
fem.nodes.constraints       # node-pair and node-group records
fem.elements.constraints    # surface interpolation records

The declaration is pure intent — no node exists yet when you make it, which is exactly why it survives remeshing. At get_fem_data(...) the resolver walks your declarations against the real mesh: it finds the co-located pairs, runs the projections, computes the interpolation weights, and lands concrete records on the snapshot — node-level constraints on fem.nodes.constraints, surface-coupling records on fem.elements.constraints.

Then the third stage happens without you: the typed apeSees(fem) bridge reads those records and emits the solver commands automatically — equalDOF, rigidLink, rigidDiaphragm, the ASDEmbeddedNodeElement penalty elements, phantom nodes and all. You declare the tie; the bridge writes the deck. The one refinement worth knowing at concept level: in a staged analysis, a constraint whose nodes only come alive in a later stage can be given a name= at declaration and claimed inside that stage's block (s.tie(name=...), s.embedded(name=...)) so it emits there instead of globally.

One convention runs through all of it: DOFs are 1-based indices (1=ux, 2=uy, 3=uz, 4=rx, 5=ry, 6=rz), and every geometric tolerance is in model units — 1e-6 is right for a metre model and uselessly tight for a millimetre one. The single exception to the index convention is g.constraints.bc(...), the fix-to-ground boundary condition, which takes an OpenSees-style restraint mask (dofs=[1, 1, 0] means "fix x and y") because it becomes ops.fix downstream. It lives on g.constraints because fixity is permanent, not load-pattern-scoped; the recipe is in Supports & BCs.

Same mesh: co-located pairs and rigid clusters

When two parts genuinely share node positions at their boundary — a conformal interface, typically the product of g.parts.fragment_all() or careful construction — the constraint is just bookkeeping. equal_dof finds every master/slave node pair within tolerance and ties the selected DOFs:

g.constraints.equal_dof(
    "slab", "column_top",
    dofs=[1, 2, 3],        # couple translations only
    tolerance=1e-6,
)

Omit dofs and all available DOFs are tied. Two siblings cover the less-common shapes of the same situation: equal_dof_mixed ties differently numbered DOFs across a pair (master uz driving slave rz, say, at a solid-to-shell joint), and penalty replaces the algebraic constraint with a stiff spring when the solver's constraint handler struggles.

rigid_link adds kinematics: each slave follows the master through a rigid offset arm, so a rotation at the master produces translations at the slaves. link_type="beam" carries all six DOFs; "rod" couples translations only. Three group-level constraints build on the same idea. rigid_diaphragm is the classic floor constraint — every slab node within plane_tolerance of the diaphragm plane shares in-plane motion with a master node:

g.constraints.rigid_diaphragm(
    "slab", "slab_master",
    master_point=(2.5, 2.5, 3.0),
    plane_normal=(0, 0, 1),
    plane_tolerance=0.05,
)

rigid_body welds a whole region to a master (all six DOFs), and kinematic_coupling is the RBE2 of the family: a reference node rigidly drives a node set with a per-set DOF selection, carrying the correct moment-arm transport so an offset reference couples rigidly. Its flexible counterpart is distributing_coupling — RBE3 — which distributes a force at a reference point over a node set as a statically equivalent pattern while the set stays deformable, with weighting="area" for a traction-like spread. Reach for RBE2 when the region must move rigidly (a loading platen); RBE3 when you want to introduce a load without stiffening anything. Both emit dedicated Ladruno fork elements: the deck is written on any build, but running it needs the fork. A labelled g.decouple_node (or its handle) is a valid reference for either — pass it as master_label so essential (ops.fix) or natural (p.load) BCs land on the ndf=6 work point without a post-mesh retarget.

Non-matching meshes: ties

Fragmenting everything conformal is not always possible or desirable — a shell floor on a solid column, two independently meshed blocks, a refined region against a coarse one. tie handles displacement continuity across a non-matching interface: each slave node is projected onto the closest master face and its DOFs are interpolated from that face's corner DOFs through the face's shape functions, u_slave = Σ Nᵢ(ξ, η) · u_masterᵢ — the same idea as Abaqus *TIE.

g.constraints.tie("flange_surface", "web_surface", tolerance=1.0)

The tolerance is the maximum projection distance — generous enough to bridge any geometric gap, tight enough not to grab the wrong face. As a rule, make the finer mesh the master: more shape functions to project onto. When a part has many faces, scope the search with master_entities= / slave_entities=.

How the tie is enforced is a separate choice, exposed as enforce=. The default "penalty" emits ASDEmbeddedNodeElement penalty elements — robust and handler-independent, with a stiffness that defaults to "auto" (sized from the host material at emit, so it lands a few orders above the element stiffness in whatever unit system you model in). "equation" emits exact multi-point equations (translations only) enforced by a constraint handler, and "penalty_al" uses the fork's augmented-Lagrange penalty element. Start with the default; switch to "equation" when you need the interface exact and you are running on a Ladruno fork build — the live equation route is fork-gated, though deck emission is not (Backend capabilities).

How the tie's weights are computed is a third, orthogonal choice — method=. The default "collocation" is the projection described above: each slave node pinned to the master field at one point. That is the right answer until the two sides differ in element order — hex20 faces tied onto hex8 faces — where pinning quadratic nodes to a bilinear field over-constrains the finer side. method="mortar" instead integrates the bond over the facet overlaps with a dual (biorthogonal) slave basis, so neither side's interpolation is imposed on the other, and the result is master/slave symmetric in a way collocation is not. It requires enforce="equation", a flat coincident interface, and dim-2 element groups on both sides (extract with dim=None); every degenerate case is a hard error rather than a silent zero-record tie. It works on composed assemblies, which is where order-mismatched models necessarily live. The practical rules are in Tie non-matching meshes.

tied_contact is the surface-to-surface version of the same projection — every slave-surface node tied to the master surface. It is one-directional (slave conforms to master; pick the finer mesh as master), and takes the same enforce= routes as tie.

The step-by-step recipe is in Tie non-matching meshes, and a complete worked model in the tie example.

Buried parts: embedded and beam-to-solid

Two situations look like ties but aren't surface-to-surface. The first is embedding: a lower-dimensional part living inside a host — rebar curves in a concrete volume, stiffeners in a shell.

g.constraints.embedded("concrete_volume", "rebar_lines")

Each embedded node is located inside a host element and constrained to follow the host's displacement field through its shape functions. Non-simplex and higher-order hosts are decomposed to linear sub-tets/sub-tris under the hood, and embedded nodes that already coincide with a host corner are dropped — they're attached through shared connectivity and constraining them again would be redundant. Emission is ASDEmbeddedNodeElement, automatic like everything else.

The second is the DOF-mismatch bridge: connecting a 6-DOF frame node to a 3-DOF solid face. node_to_surface builds a compound constraint — phantom 6-DOF nodes duplicated at the slave positions, rigid links from the master to each phantom, equalDOF from each phantom down to the real solid node:

g.constraints.node_to_surface("frame_end", "solid_face")

Here the master is a geometric point entity, not a part label — the one place in the family where you target an entity rather than a name. The node_to_surface_spring variant replaces the rigid links with stiff beam elements; use it when the master carries free rotational DOFs under direct moment loading, where the purely kinematic version leaves those rotations without a stiffness path and the matrix conditioning suffers.

Contact

Everything above is kinematic — a permanent bond, active from the first step. g.constraints.contact(...) is the genuinely different animal: a face-to-face contact interaction that can open, close, slide, and carry friction, emitted through the Ladruno fork's contact subsystem.

g.constraints.contact("wall_face", "soil_face",
                      formulation="nts", kn="auto", mu=0.4)

Two formulations: "nts" (node-to-segment penalty, with Coulomb friction) and "mortar" (segment-to-segment augmented-Lagrange — the accuracy lane for non-matching interfaces). Passing tie=True to the mortar formulation freezes the pair into a permanent contact-tie bond. For a plain permanent bond, prefer tie(method="mortar", enforce="equation") instead: the contact-tie is enforced by penalty/augmentation (the gap is driven toward zero, never exactly zero), it mandates the LadrunoContact handler — which cannot coexist with any enforce="equation" tie in the same analysis — and it accepts linear facets only. tie=True earns its keep when the interface must also behave as contact later in the run. Contact declarations resolve additively onto fem.elements.contacts rather than the MP-constraint channels; they need a live gmsh session (declaring one on a from_h5/composed session raises), and like the fork coupling elements, the deck emits on any build but runs only on the fork.

One historical note, since older material mentions it: g.constraints.mortar() still exists but is a deprecated alias that delegates to contact(formulation="mortar", tie=True) and warns. Call contact() directly.

One more lane has no master mesh at all. g.constraints.contact_plane puts the slave surface against a fixed infinite rigid plane given by a normal and a point, frictionless, with a kn you must supply — there is no "auto" here. It is the cheapest real contact in the library, and the only one that does not demand ndf == ndm, so a shell can sit on a floor.

Contact in a plane model

Everything above holds in 2D, but the shape of a contact surface changes and so does the way you orient it. The fork decides a surface's dimension by reading its nodes' coordinates rather than the interpreter's ndm, and apeGmsh writes exactly ndm coordinates per node, so a plane model reaches the 2D lanes without being told. What changes is what you name: a 2D contact surface is the meshed boundary curve — a dim-1 physical group — and not the dim-2 body. Naming the body is refused rather than quietly turned into a contact between solid elements.

Underneath, the fork's 2D surface is a flat list of node pairs chained head-to-tail, so three segments need six tags rather than four. The four-tag shorthand is legal to the parser and declares two disjoint segments with a hole between them — a deck that converges, balances its reactions, and transmits the load through the wrong distribution. apeGmsh builds that chain from the physical group's own edge connectivity, on both sides, so the holed form cannot be produced. That is the main reason to route a plane contact through the library rather than write the deck by hand.

Orientation is where 2D genuinely differs. In 3D the kernel derives a correct normal per facet, which is why apeGmsh never guesses a global one there. The 2D lanes instead take a single interface-level sign from a centroid vote, and that vote is ambiguous the moment the two surfaces are coincident — which in 2D is the ordinary case rather than the exception: the masonry joint, the footing seated on soil, every zero-gap interface. A flush declaration is refused by name, and there are two ways to answer it. outward=(ox, oy) names a direction toward the slave's allowed half-space and works on both formulations. outward="winding" instead declares the side through the master chain's own travel — the slave lies to the left of it — which is the only thing that can orient a curved or closed master, since no single direction can. Winding is available on the NTS formulation only: the mortar lane runs no chain-integrity scan for it to rest on, so a flush mortar interface always takes the vector and a curved mortar master cannot be declared at all. What neither side will do is guess — a master surface facing away from its slave is refused too, because the fork's vote picks a sign and never asks whether this was the face you meant.

The remaining 2D-only parameter is thickness, worth stating precisely because a plane model carries an out-of-plane thickness in three unrelated places. The element's own thickness is baked into element stiffness and contact never re-reads it. The mortar thickness=h scales the explicit penalties and the tie stiffness once, at the fork's injection site; it is mortar-only, refused on NTS and in 3D. And eps_n="auto" is deliberately left alone — it already absorbs the element thickness through the element's initial stiffness, so scaling it again would be a squared error.

Three properties of the model decide whether a plane contact transmits anything at all, and none of them is visible in the declaration. A zero initial gap does not arm the NTS lane, so seed a small overlap — the rigid plane arms from zero and needs none. 2D contact carries the normal direction only, so a body whose only other support is the interface keeps a free transverse mode, and when that mode is excited Newton drifts while the solver reports a misleading warning about the interface geometry. And the facets of a curved or closed master must be sized from the expected penetration, or from the loop's own extent, rather than from the elastic mesh around them: drive them from the same size parameter and refining the model makes the contact worse, while a coarse closed ring converges, balances, and transmits exactly zero. Parallel 2D contact does not exist — it is out of scope in the fork, and refused by name here the moment the model is partitioned. Kernel behaviour beyond all this — the vertex policy, the radial end-cap, the units table — belongs to the fork's own LadrunoContact2D_guide.md.

Unilateral interface springs

Between a bond and full contact sits a third thing, and the model that paid for it is a tunnel liner in squeezing rock. g.constraints.interface(...) puts one spring per coincident node pair across a 2D continuum boundary — unilateral in the normal direction, capped at a bond strength in the tangential one. It needs no contact search and no exclusive handler, because the pairing is plain co-location, the same matching equal_dof does. What separates it from equal_dof, tie and embedded is that it is not a bond at all: the normal law carries compression and nothing else, and the tangential law stops carrying shear once τ_b × A_trib is reached.

That distinction is mechanical, not cosmetic. With a bilateral bond the model cannot reproduce demand saturation — the ground converges, the tie transmits whatever force that convergence implies, and the liner's demand grows without bound, while the field record shows metre-class convergence coexisting with damaged-but-standing arches. Capping the bond puts a ceiling on what the ground can hand the structure, and letting the interface open means it carries nothing where the two sides have parted. The surface, including the sign convention you should check once per model, is in the constraints API; the two things to get right at declaration are the out-of-plane thickness (required, never guessed) and slave_ndf= when the coincident wire will be meshed as a beam.

Reading back what you declared

After resolution the snapshot can tell you what actually happened — fem.inspect.constraint_summary() prints the record counts by kind, and fem.nodes.constraints.summary() / fem.elements.constraints.summary() give DataFrame views. If a tie found zero pairs, the answer is almost always the tolerance: too tight finds nothing, too loose couples nodes that shouldn't be coupled. Check the summary before blaming the solver.

---

Next: Loads & masses.