Design Note: Rigid Element Formulation for node_to_surface¶
Status: Aspirational — design note only; the in-use workaround is
node_to_surface_spring(stiff-beam variant).
Problem¶
rigidLink in OpenSees is an MP_Constraint (constraint equation), not an element.
When the Transformation handler eliminates slave DOFs, it substitutes
u_slave = u_master + theta_master x r into the global equations. But this
never adds stiffness to the master's rotation DOFs. If no element connects
directly to those DOFs, the diagonal of K is zero there — the system is singular.
This manifests as: - Newton failure under moment loading at the reference node - Singular K even with no load if rotation DOFs are free and unconnected
The current workaround (node_to_surface_spring) replaces rigidLink with a
stiff elasticBeamColumn. This works because the beam's 12x12 K matrix gives
direct stiffness to both translation and rotation DOFs at both ends. But it is
approximate — the beam has finite (though very large) stiffness.
Proposed Solution: Rigid Body Element¶
Formulate the rigid-body kinematics as a finite element rather than a constraint equation. The element encodes:
directly in its element stiffness matrix, so K assembly gives the master rotation DOFs proper diagonal entries.
Two implementation paths¶
A. Penalty rigid element¶
An element with K proportional to a penalty parameter alpha:
where G is the constraint matrix that maps master DOFs to slave displacements. For a single master-slave pair with arm vector r = (rx, ry, rz):
The element residual is f = K_e * (u_slave - u_master - theta x r).
Pros: simple, no extra DOFs, direct stiffness on all master DOFs. Cons: penalty sensitivity (too large -> ill-conditioned, too small -> constraint violation).
B. Lagrange multiplier rigid element¶
Add multiplier DOFs lambda that enforce the constraint exactly:
Pros: exact constraint, no penalty tuning. Cons: adds DOFs, saddle-point system needs compatible solver, zero diagonal block.
C. Augmented Lagrangian¶
Combine penalty + Lagrange: K_aug = K_struct + alpha * G^T G, plus multiplier update. Pros: exact in the limit, well-conditioned. Cons: iterative multiplier update, more complex.
Recommendation¶
Path A (penalty) is the pragmatic choice for OpenSees:
- OpenSees already uses penalty for other constraint types
- No extra DOFs or solver changes
- Penalty value can be derived from material/geometry (same as stiff beam E*A/L)
- Can be wrapped as a custom element via OpenSeesPy element('genericClient', ...)
or implemented as a Python-side K matrix contribution
The stiff elasticBeamColumn in node_to_surface_spring is already a penalty
rigid element in disguise — its stiffness is the penalty. A dedicated element
would be cleaner (no need for vec_xz, no geometric length sensitivity) but
mechanically equivalent.
Scope for Implementation¶
- Custom element class in apeGmsh that computes K_e from the arm vector r and a penalty stiffness derived from the real structure's material properties
- Emit via
element()call in the OpenSees writer — either as agenericClientelement or by assembling the K contribution externally - Replace
node_to_surface_springstiff beam emission with the rigid element - Fallback: keep stiff beam as the default since it works and is well-tested; offer rigid element as an option
Relationship to Existing Code¶
node_to_surface(rigidLink): fails under free-rotation + moment loadingnode_to_surface_spring(stiff beam): works, approximate, vec_xz sensitivitynode_to_surface_rigid_element(proposed): works, exact or near-exact, no vec_xz
The phantom node topology remains the same in all three variants:
Only the [link] implementation changes.
When the Issue Does NOT Apply¶
- SP-prescribed rotation DOFs: eliminated from system, no singularity
- Master connected to frame elements: frame beam gives rotation stiffness
- Pure force on translations only: rotation DOFs are zero but may still be singular (no stiffness); works in practice only if the solver tolerates the zero pivot (some do, some don't)