Element-kernel authoring#

How to add a new element kernel — shape functions, stiffness / mass assembly, real-constant layout, registration. Companion to Interop reader authoring (the reader-side of “a deck mentions CBUSH, our parser doesn’t know it”). This page is the kernel-side: “the parser can route to BUSH6, but BUSH6 doesn’t exist yet.”

When you’re authoring a new kernel #

Two driving paths get you here.

A VM hits a card whose target element type isn’t in the kernel registry. The reader gap is filed under the [interop] label, the kernel gap under [kernel]; both carry verify-blocked until they land. See the Coordinated card+kernel landings section of Interop reader authoring for the joint workflow.
A new element formulation (EAS, MITC, hourglass-stabilised reduced integration, …) is needed to unlock a benchmark family that the existing kernel can’t reach. Author as a sibling module under src/femorph_solver/elements/ and register under a distinct neutral name; route the foreign-deck KEYOPT / SECTYPE to it via the spec class.

Whichever path got you here, the kernel itself is independent of how the reader routes to it. The contract below is what makes a kernel pluggable.

The contract — what a kernel must provide #

Every kernel subclasses femorph_solver.elements._base.ElementBase and provides:

@register                                      # registers under cls.name
class Spring(ElementBase):
    name = "SPRING"                            # neutral kernel name
    n_nodes = 2                                # connectivity arity
    n_dof_per_node = 3                         # 3 = trans-only, 6 = trans+rot
    vtk_cell_type = 3                          # VTK_LINE / VTK_QUAD / VTK_HEXAHEDRON / ...
    # dof_labels defaults to CANONICAL_DOFS[:n_dof_per_node]

    @staticmethod
    def ke(coords, material, real):            # element stiffness, dense (n_dof, n_dof)
        ...
    @staticmethod
    def me(coords, material, real, lumped=False):  # element mass
        ...

The base class fills dof_labels from CANONICAL_DOFS = ("UX", "UY", "UZ", "ROTX", "ROTY", "ROTZ")[:n_dof_per_node] unless the subclass overrides it (e.g. shells with explicit drilling DOFs).

The neutral-name convention #

Kernels register under their topology-first neutral name:

Neutral	MAPDL alias	NASTRAN	Abaqus
`HEX8`	`SOLID185` / `SOLID45`	`CHEXA` (8 nodes)	`C3D8` / `C3D8I` / `C3D8R`
`HEX20`	`SOLID186` / `SOLID95`	`CHEXA` (20 nodes)	`C3D20`
`TET10`	`SOLID187` / `SOLID92`	`CTETRA` (10 nodes)	`C3D10`
`QUAD4_SHELL`	`SHELL181` / `SHELL63`	`CQUAD4` w/ `PSHELL`	`S4` / `S4R`
`QUAD4_PLANE`	`PLANE182`	`CQUAD4` w/ `PSHELL MID2=0`	`CPS4` / `CPE4`
`BEAM2`	`BEAM188`	`CBAR` / `CBEAM`	`B31`
`TRUSS2`	`LINK180` / `LINK8`	`CROD` / `CONROD`	`T3D2`
`POINT_MASS`	`MASS21`	`CONM2`	`*MASS`

Foreign-deck name → neutral name happens at the interop boundary (see Interop reader authoring’s _ELEMENT_MAP / _ELEMENT_TYPE_MAP). The kernel never sees a vendor name.

Real-constant layout #

The real tuple’s slot meaning is per-element-family and the reader’s responsibility. Document the layout at the top of the kernel module so future readers can audit cross-vendor parity:

"""Spring — 3D 2-node spring-damper (longitudinal mode only).

Real constants:

    real[0] : K    — spring stiffness (force / length)
    real[1] : CV1  — linear damping (force · time / length, unused for K_e)
    real[2] : CV2  — non-linear (cubic) damping (unused here)
    real[3] : IL   — initial length (unused)
"""

For beam-family elements specifically, the layout is fixed at (A, IZZ, IYY, J) — see Interop reader authoring “The materialisation contract” for the reader-side I1/I2 swap that gets you there.

Material dictionary #

material is a dict keyed by canonical property names:

EX — Young’s modulus.
PRXY — Poisson’s ratio.
DENS — density.
ALPX — thermal expansion coefficient.
GXY (optional) — shear modulus override; default E / (2(1 + ν)).

Per-formulation hints (e.g. plane stress vs plane strain, EAS vs B-bar) flow in via _QUAD4_PLANE_MODE, _QUAD4_PLANE_TECH underscore-prefixed keys — the spec instance writes them.

DOF ordering #

The element returns ke and me in node-major, DOF-minor order:

[u₁ⁿ⁰, u₂ⁿ⁰, …, u_{n_dof}ⁿ⁰,   u₁ⁿ¹, u₂ⁿ¹, …, u_{n_dof}ⁿ¹,   …]

i.e. the first n_dof_per_node rows are node 0’s DOFs in dof_labels order, then node 1’s, etc. The assembly path relies on this; do not ship a kernel with arbitrary ordering.

The sequence — adding a new kernel #

Worked example: a hypothetical BUSH6 element — 2-node 6-DOF spring/damper connector with three translational and three rotational stiffnesses.

Step 1 — write the kernel module #

Drop src/femorph_solver/elements/bush6.py. Header docstring cites the published source for the formulation (this is Provenance inventory’s job — every non-trivial numerical kernel must trace to a paper / textbook).

"""BUSH6 — 2-node 6-DOF translational+rotational connector.

Real constants:

    real[0..5]  : K_t1, K_t2, K_t3, K_r1, K_r2, K_r3
                  (per-DOF stiffnesses in element local frame)

References
----------
* Cook, Malkus, Plesha, Witt, *Concepts and Applications of
  Finite Element Analysis*, 4th ed., Wiley, 2002, §2.10
  (multi-DOF spring-damper between two nodes).
"""

from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register


@register
class Bush6(ElementBase):
    name = "BUSH6"
    n_nodes = 2
    n_dof_per_node = 6
    vtk_cell_type = 3  # VTK_LINE

    @staticmethod
    def ke(coords, material, real):
        r = np.asarray(real, dtype=np.float64)
        if r.size < 6:
            raise ValueError("Bush6 requires REAL[0..5] = (K_t1, K_t2, K_t3, K_r1, K_r2, K_r3)")
        # Local stiffness — diagonal 6×6 then sandwiched between two nodes.
        k_local = np.diag(r[:6])
        # Block layout: [[ K, -K], [-K, K]] in 12×12 element frame.
        K = np.block([[k_local, -k_local], [-k_local, k_local]])
        return np.ascontiguousarray(K)

    @staticmethod
    def me(coords, material, real, lumped=False):
        return np.zeros((12, 12), dtype=np.float64)

Conventions:

ke and me are static methods. Kernels carry no instance state; everything they need comes through the three arguments.
Validation lives at the kernel boundary — raise with a clear message naming the missing field. The interop layer should have already populated real correctly, but the kernel shouldn’t trust it.
Return np.ascontiguousarray(...) for stiffness so the assembly path doesn’t pay for a copy.

Step 2 — register the spec #

Append to src/femorph_solver/elements/_specs.py:

@dataclass(frozen=True)
class Bush6Spec(ElementSpec):
    name: str = "BUSH6"
    # Add formulation kwargs here if the element has them
    # (e.g. tech="enhanced" on QUAD4_PLANE).

This lets Model.assign(Bush6Spec, ...) and the spec-by-name lookup work. For the simplest cases (no formulation kwargs) the spec is just the name + class binding.

Step 3 — wire the importer #

Add to src/femorph_solver/elements/__init__.py so the kernel is registered at import time:

from femorph_solver.elements.bush6 import Bush6  # noqa: F401

Without the import, @register never runs and elements.get("BUSH6") raises.

Step 4 — element-kernel unit tests #

Drop tests/elements/bush6/test_bush6_kernel.py with the element-only invariants:

Stiffness symmetry: np.allclose(K, K.T).
Rigid-body modes: K @ rbm == 0 for the 6 rigid-body modes expected for the topology.
Patch-test or single-element analytical solution: a 2-node bush under prescribed displacement returns the analytical reaction.
Real-constant layout: each slot routes to the right physical effect.

Don’t run the assembled solver here — that’s the cross-solver harness’s job. Element tests live in tests/elements/<element>/ per Test-suite layout.

Family	Reference module	Idiomatic shape
3D solid	`hex8.py` / `tet10.py` / `hex20.py`	Gauss-quadrature loop over the canonical reference element; `B = ∂N/∂x` via Jacobian; `K += Bᵀ D B det(J) w`. EAS / hourglass control add an extra enrichment step inside the same loop.
2D plane / shell membrane	`quad4_plane.py` / `quad8_plane.py`	2×2 Gauss for membrane; thickness multiplies the constitutive `D`; plane stress vs strain selects which 3×3 `D` is built.
Shell (Mindlin-Reissner)	`quad4_shell.py`	Membrane + bending + transverse shear. Reduced integration on the shear term (1×1 Gauss) is standard; drilling stabilisation (Hughes-Brezzi penalty) keeps ROTZ from being a free DOF.
Beam (Hermite cubic)	`beam2.py`	Local frame from orientation vector → 12×12 local K / M → transform to 6-DOF-per-node global frame via block rotation matrix.
Rod / truss	`truss2.py`	1-D constant-strain element; `K = (EA/L) (d ⊗ d)` sandwiched between two nodes.
Spring / point connector	`spring.py`	Same direction-cosine sandwich pattern as truss with a scalar real instead of `EA`.
Lumped mass	`point_mass.py`	`ke` is zero; `me` is diagonal with the real-constant mass on translational DOFs and rotary inertia on rotational DOFs.

Common pitfalls #

Forgetting the import in elements/__init__.py. @register only runs when the module is imported. A kernel that’s never imported is invisible to the registry.
Returning a non-symmetric ke. The assembly path adds ke directly into a CSR K; non-symmetry surfaces only at modal-solve time as complex eigenvalues. Test for symmetry in the kernel unit test, every time.
Wrong DOF ordering. The “node-major, DOF-minor” rule isn’t enforced by the type system; eyeballing the 12×12 (or 24×24) matrix is the only safety net. Cross-check against a 1-element analytical solution before trusting it.
Hard-coded shape-function constants instead of building them from the canonical reference element. hex20 has three families of nodes (corner / mid-edge / mid-face) with different shape-function forms — copy the existing kernel’s pattern, don’t re-derive.
Mutating coords / material / real. Kernels are pure functions; the assembly loop reuses these arrays across many elements. In-place edits cascade silently.
Mass matrix formulation. Default is consistent mass; lumped=True switches to row-sum lumping. Modal solves (Model.solve_modal) use consistent unless the problem explicitly requests lumped — the closed-form references usually assume consistent.
Missing the References section in the docstring. Every non-trivial kernel must trace to a paper / textbook. See Provenance inventory; audit rows in that table are flagged for follow-up.

Where things live #

Concern	Path
Element kernel module	`src/femorph_solver/elements/<name>.py`
Element registry	`src/femorph_solver/elements/_registry.py`
Element specs (formulation kwargs, KEYOPT routing)	`src/femorph_solver/elements/_specs.py`
Base class	`src/femorph_solver/elements/_base.py`
Per-element unit tests	`tests/elements/<element>/test_<element>_*.py`
Beam-section helpers	`src/femorph_solver/elements/_beam_sections.py`
Stress-recovery helpers	`src/femorph_solver/elements/_stress.py` and `src/femorph_solver/result/_stress_recovery.py`
Provenance entry	`doc/source/dev/provenance.rst` (file’s `References:` block + the matrix row)

Companion pages: Interop reader authoring (reader-side of the “new card → new kernel” pair), Verification-manual spec (landing the kernel against a published benchmark), Provenance inventory (kernel-citation inventory).