Element-kernel authoring ======================== How to add a new element kernel — shape functions, stiffness / mass assembly, real-constant layout, registration. Companion to :doc:`interop_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." .. contents:: Page contents :local: :depth: 2 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 :doc:`interop_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 :class:`femorph_solver.elements._base.ElementBase` and provides: .. code-block:: python @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**: .. list-table:: :header-rows: 1 :widths: 22 26 26 26 * - 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 :doc:`interop_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: .. code-block:: python """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 :doc:`interop_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 :doc:`provenance`'s job — every non-trivial numerical kernel must trace to a paper / textbook). .. code-block:: python """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``: .. code-block:: python @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: .. code-block:: python 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//`` per :doc:`testing`. Step 5 — extend the interop maps ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If the new kernel has a vendor analogue, route the foreign-deck name to it. See :doc:`interop_authoring`'s "Adding a card to ``from_bdf``" Step 4 for the ``_ELEMENT_MAP`` / ``_ELEMENT_TYPE_MAP`` edits. Step 6 — wire a verification row ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The kernel earns its keep against a published benchmark. Land the registry row + fixture pair in the same PR (or as a gating-on-each-other follow-up) — see :doc:`verification_manual_spec` for the registry entry shape. Patterns by element family -------------------------- The reference kernels in ``src/femorph_solver/elements/`` cover the standard families. When you author a new one, mirror the closest existing kernel's idioms — the pattern below should land: .. list-table:: :header-rows: 1 :widths: 22 26 52 * - 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 :doc:`provenance`; ``audit`` rows in that table are flagged for follow-up. Where things live ----------------- .. list-table:: :header-rows: 1 :widths: 32 68 * - Concern - Path * - Element kernel module - ``src/femorph_solver/elements/.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//test__*.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: :doc:`interop_authoring` (reader-side of the "new card → new kernel" pair), :doc:`verification_manual_spec` (landing the kernel against a published benchmark), :doc:`provenance` (kernel-citation inventory).