FESADev/docs/linear-static-3d-euler-beam/numerical-review.md

# Linear Static 3D Euler Beam Numerical Review

## Metadata

- feature_id: `linear-static-3d-euler-beam`
- source_formulation: `docs/linear-static-3d-euler-beam/formulation.md`
- source_requirements: `docs/linear-static-3d-euler-beam/requirements.md`
- source_research: `docs/linear-static-3d-euler-beam/research.md`
- approved_design: `docs/superpowers/specs/2026-08-08-linear-static-3d-euler-beam-design.md`
- status: `pass-for-implementation-planning`
- owner_agent: `numerical-review-agent`
- date: `2026-08-09`
- scope: formulation correctness and implementation-planning readiness only

## Review Verdict

- verdict: `pass-for-implementation-planning`
- reason: The DOF/sign convention, natural boundary actions, 4x12 B matrix, 12x12 stiffness, 2-point Gauss rule, transformation, constant line-load vector, constrained partition, and residual reaction are mutually consistent. No confirmed mathematical defect was found.
- qualification: This verdict permits implementation planning only; it is not reference-comparison, physics-sanity, release, or full Abaqus-compatibility approval.

## Critical Findings

### Confirmed defects

- None.

### Independent DOF, curvature, and end-action derivation

The fixed local nodal order is

$$
[u_1,v_1,w_1,\theta_{x1},\theta_{y1},\theta_{z1},
 u_2,v_2,w_2,\theta_{x2},\theta_{y2},\theta_{z2}].
$$

For a small right-handed section rotation,
$\boldsymbol\theta\times[0,y,z]^T
=[z\theta_y-y\theta_z,-z\theta_x,y\theta_x]^T$.
The zero-shear constraints are therefore
$\theta_z=v'$ and $\theta_y=-w'$. Hence

$$
\varepsilon_{xx}=u'-zw''-yv'
=\varepsilon_0+z\kappa_y-y\kappa_z,
\quad \kappa_y=-w'',\quad \kappa_z=v''.
$$

Using the reviewed section-force definitions independently gives

$$
M_y=EI_y\kappa_y=-EI_yw'',\qquad
M_z=EI_z\kappa_z=EI_zv''.
$$

Twice integrating the bending virtual work, with outward sign $n=-1$ at
$x=0$ and $n=+1$ at $x=L$, yields

$$
F_y=-nM_z',\quad M_z^{end}=nM_z,\qquad
F_z=nM_y',\quad M_y^{end}=nM_y.
$$

The axial and torsional pairs are $F_x=nN$ and $M_x=nT$. Direct expansion of
$K_ld_l$ reproduces these relations at both endpoints when the V0 distributed
load is zero. Thus the equilibrium end action is an outward action, while
$DBd$ is the positive-local-x section resultant before the outward-normal
sign is applied. These outputs must remain distinct.

### B matrix, stiffness, quadrature, and modes

- Differentiating the Hermite fields reproduces the documented 4x12 B matrix. The $\kappa_y=-w''$ row has $[-H_1'',H_2'',-H_3'',H_4'']$ in the $[w_1,\theta_{y1},w_2,\theta_{y2}]$ columns; the $\kappa_z=v''$ row has $[H_1'',H_2'',H_3'',H_4'']$ in the $[v_1,\theta_{z1},v_2,\theta_{z2}]$ columns.
- Axial/torsional B rows are constant and bending rows are linear in $\xi$, so $B^TDBJ_x$ is degree at most two. The 2-point Gauss rule, exact through degree three, exactly produces the closed-form 12x12 stiffness.
- The axial and torsion blocks each have rank one, and each exact bending block has rank two. With positive $EA,GJ,EI_y,EI_z$, the free-element stiffness is symmetric positive semidefinite with rank 6.
- The six rigid modes are the three equal nodal translations, equal $\theta_x$, $\theta_{y1}=\theta_{y2}$ with $w_2-w_1=-L\theta_y$, and $\theta_{z1}=\theta_{z2}$ with $v_2-v_1=L\theta_z$. Each gives $Bd=0$.
- Since $d^TKd=\int(Bd)^TD(Bd)dx$, energy is strictly positive outside the six-dimensional rigid kernel, but the unconstrained element itself is not positive definite.
- A 1-point bending rule reduces the total element rank to 4 and introduces two spurious zero-energy modes. Production must use the approved 2-point rule or the reviewed closed form.

### Transformation and energy invariance

Projection/normalization of `n1` gives orthogonal unit $e_x,e_y$, and
$e_z=e_x\times e_y$ makes $R=[e_x^T;e_y^T;e_z^T]$ orthogonal and
right-handed: $RR^T=I$, $\det R=+1$. With
$T=\operatorname{diag}(R,R,R,R)$ in node translation/rotation block order,

$$
d_l=Td_g,\quad K_g=T^TK_lT,\quad f_g=T^Tf_l,
\quad d_g^TK_gd_g=d_l^TK_ld_l.
$$

This also confirms virtual-work invariance. Reversing the transform direction
would violate the documented component contract.

### Constant local line-load vector

Exact integration of $N_q^T[p_x,p_y,p_z,m_x,0,0]^T$ gives equal half-span
nodal forces/torques and

$$
M_{y1}=-p_zL^2/12,\quad M_{y2}=+p_zL^2/12,\qquad
M_{z1}=+p_yL^2/12,\quad M_{z2}=-p_yL^2/12.
$$

The signs follow from the $-H_2,-H_4$ interpolation for $w$ and
$H_2,H_4$ for $v$. The 2-point rule integrates the cubic load integrand
exactly. This remains a formulation-only kernel: `*DLOAD` input, a Domain
distributed-load object, and CLI distributed loading are out of scope.

### Independent numerical cross-check

An independent double-precision evaluation at $L=3.7$ with unequal positive
$EA,GJ,EI_y,EI_z$ produced:

| check | observed result |
| --- | ---: |
| 2-point Gauss/closed-form normalized error | $3.20\times10^{-16}$ |
| symmetry normalized error | $2.00\times10^{-17}$ |
| six rigid-mode residual norms | 0 to $2.81\times10^{-15}$ |
| numerical rank at relative $10^{-10}$ | rank 6 |
| constant line-load normalized error | $3.08\times10^{-16}$ |
| 1-point stiffness | rank 4 |
| rotated local/global energy difference | 0 |
| $\det R$ | $0.9999999999999999$ |

These calculations corroborate the derivation but do not replace production
tests.

## Numerical Risk Assessment

Confirmed defects, risks, and open issues are separated.

| id | risk or limit | assessment and required control |
| --- | --- | --- |
| NR-R01 | zero/near-zero length | Singular Jacobian and divergent $L^{-3}$ terms; `NR-T05` tests the approved scale-aware boundary before evaluation. |
| NR-R02 | zero/tangent-parallel `n1` | Undefined axes and inertia directions; `NR-T05` tests the projection boundary and forbids a silent fallback. |
| NR-R03 | under-integration | 1-point integration adds two mechanisms; `NR-T03` must show 2-point rank 6 and negative-control rank 4. |
| NR-R04 | mixed translation/rotation scaling | Raw rank is length-unit dependent; `NR-T04` uses length scaling and a well-scaled fixture. |
| NR-R05 | incomplete constraints | `Kff` is singular while any rigid mode survives; `NR-T09` covers free, partial, and stable systems. |
| NR-R06 | extreme stiffness ratios/units | A mathematically SPD system may be ill-conditioned; `NR-T10` checks residual/analytical error. Automatic scaling or a new threshold is not approved. |
| NR-R07 | transform direction | Rotated bending signs/components can be wrong; `NR-T05` checks orthogonality, handedness, work, and energy. |
| NR-R08 | end-action/section-cut confusion | Opposite signs can contaminate comparison; `NR-T07` checks both endpoints and distinct result types. |
| NR-R09 | sparse duplicate reduction | Thread-dependent summation can break reproducibility; `NR-T11` checks canonical CSR and normalized values. |
| NR-R10 | short/deep beam | Euler–Bernoulli may be too stiff; explicit limitation only. No unapproved slenderness diagnostic and no B31 reinterpretation. |
| NR-R11 | warping, eccentricity, or `I12!=0` | The diagonal constitutive model is insufficient; explicit rejection/exclusion, never silent approximation. |
| NR-R12 | result averaging | Averaging can hide discontinuities; baseline nodal averaging is excluded and `NR-T07` preserves endpoint identity. |

- shear_locking: N/A because no independent shear strain exists; deep-beam applicability remains limited.
- volumetric_locking: N/A.
- hourglass: N/A for the approved 2-point rule; the distinct 1-point rank deficiency is NR-R03.
- distortion: no curved/distorted mapping exists in this straight affine element; length and axis singularities are NR-R01/NR-R02.
- nonlinear convergence: N/A; linear residual and conditioning are NR-T09/NR-T10.
- mesh convergence: tip-force/tip-moment polynomial fields are one-element exact; the formulation-only constant transverse load has a quartic exact displacement and requires refinement evidence.

## Consistency Checks

| check | result | evidence |
| --- | --- | --- |
| units | pass | Stiffness blocks and nodal work pairs are dimensionally consistent. |
| dimensions | pass | B is 4x12, D is 4x4, and $B^TDBJ_x$ is 12x12. |
| signs | pass | Cross-product, curvature, natural-boundary, end-action, and line-load derivations agree. |
| dof_ordering | pass | Node blocks preserve `[u,v,w,theta_x,theta_y,theta_z]`. |
| coordinate_transforms | pass | R/T are right-handed orthogonal maps with local = transform times global. |
| jacobian | pass with validation | $J_x=L/2>0$ for accepted elements. |
| constitutive | pass | Positive diagonal $EA,GJ,EI_y,EI_z$ gives nonnegative energy; `I12!=0` is rejected. |
| integration | pass | The 2-point locations, unit weights, and $J_x$ integrate stiffness and constant-load kernels exactly. |
| symmetry | pass | Variational and numerical checks agree. |
| rank 6 and rigid modes | pass | Independent ranks $1+1+2+2=6$ and six rigid vectors span the kernel. |
| positive deformation energy | pass | Strictly positive outside the rigid kernel. |
| constrained partition | pass | $K_{ff}d_f=F_f-K_{fc}d_c$ includes nonzero $d_c$. |
| residual reaction | pass | $R_c=(Kd-F)_c=K_{cf}d_f+K_{cc}d_c-F_c$. |
| output locations | pass with downstream contract | Gauss values, section resultants, end actions, and residual reactions remain distinct. |

## Verification Readiness

| test id | concrete criterion |
| --- | --- |
| `NR-T01-DOF-SIGN` | Compare `B*d` with independently differentiated Hermite fields at multiple $\xi$; normalized error $\le10^{-12}$ and explicit $\theta_y=-w'$, $M_y=-EI_yw''$. |
| `NR-T02-PATCH` | Check endpoint value/slope, constant axial strain/twist, both constant-curvature patches, and six rigid modes; rigid residual $\le10^{-10}$. |
| `NR-T03-GAUSS-CLOSED` | All 12x12 entries and symmetry normalized $\le10^{-12}$; negative 1-point control must be rank 4. |
| `NR-T04-RANK-ENERGY` | Use $Q=\operatorname{diag}(1,1,1,L,L,L,1,1,1,L,L,L)$ and $\widehat K=Q^{-T}KQ^{-1}$. For a documented well-scaled fixture require six singular values $\le10^{-10}\sigma_{max}$, six $>10^{-10}\sigma_{max}$, and positive energy for six deformation vectors. |
| `NR-T05-TRANSFORM` | For a non-axis-aligned member require normalized orthogonality, handedness, virtual-work, transformed-stiffness, and energy errors $\le10^{-12}$; test geometry thresholds. |
| `NR-T06-LINE-LOAD` | Compare all 12 signed closed-form components at normalized $\le10^{-12}$; separately verify `*DLOAD` remains unsupported. |
| `NR-T07-END-SIGNS` | Pure axial/torsion/two-plane bending states must match outward natural actions and positive-face resultants at both endpoints, normalized $\le10^{-12}$; no averaging. |
| `NR-T08-ANALYTICAL` | Axial, torsion, local-y, and local-z cantilevers at relative $\le10^{-9}$; constant transverse load requires load-vector equality and mesh convergence. |
| `NR-T09-PARTITION-SPD` | Nonzero $d_c,F_c$: displacement and free/reaction residual normalized $\le10^{-10}$. Stable `Kff` factorizes; surviving rigid modes fail structurally with no result. |
| `NR-T10-CONDITIONING` | Sweep documented valid scales; accepted well-conditioned cases require residual $\le10^{-10}$ and analytical error $\le10^{-9}$; numerically unresolved cases fail explicitly. |
| `NR-T11-DETERMINISTIC-ASSEMBLY` | Across repeated thread counts require identical CSR structure/reduction order and values normalized $\le10^{-12}$. |

### Open issues

1. `NR-O01-PARDISO-CONTRACT`: official oneMKL evidence is still required for SPD matrix type, phases, zero-based CSR, repeated RHS, conditioning/error reporting, and failure codes before that adapter plan is approved. This is non-blocking for element mathematics and connects to NR-T09/NR-T10.
2. `NR-O02-DETERMINISTIC-REDUCTION`: stable COO sort and duplicate-summation rules are project policy and must be made explicit before NR-T11.
3. `NR-O03-STATION-NORMALIZATION`: reversed connectivity/local-axis orientation and legitimate jumps at loaded interior nodes need an explicit downstream row-normalization/eligibility rule. The legacy baseline may use its documented stable orientation and unloaded interior stations, but mismatch must never be averaged. NR-T07 covers element signs.

2026-08-18 amendment: `NR-O03-STATION-NORMALIZATION` is resolved and superseded for the
approved case. The regenerated elemental-force CSV supplies `Element Label`, so comparison uses
direct `(instance, element label, endpoint node label, component)` identity and never collapses,
selects or averages adjacent endpoints.

No open issue requires formulation revision. NR-O01/NR-O02 are implementation-planning
handoffs; the resolved NR-O03 identity is fixed by the I/O and reference-model contracts.

## Required Revisions

### Formulation Agent

- None. No confirmed mathematical defect or missing derivation blocks implementation planning.

### Research Agent

- Resolve NR-O01 from official Intel documentation before finalizing the PARDISO adapter plan.

### Reference Model Agent

- Enforce the resolved NR-O03 direct element-endpoint identity without modifying reference artifacts.

## Downstream Handoff

### Implementation Planning Agent

- Convert NR-T01 through NR-T11 into TDD `RED -> GREEN -> VERIFY` work without changing approved tolerances.
- Preserve stiffness/partition/factorization before load/effective-RHS/substitution.
- Treat SPD as a post-constraint property and retain singularity/conditioning diagnostics.
- Resolve NR-O01 and NR-O02 before the corresponding implementation steps.

### I/O Definition Agent

- Preserve outward endpoint action, positive-face section resultant, Gauss result, and assembled residual reaction as distinct identities.
- Define NR-O03 without broadening the approved Abaqus subset.

### Reference Model Agent

- Cover axial, torsion, both bending planes, rotated-space transformation, endpoint signs, and nonzero prescribed displacement.
- Keep `reference/cantilever beam/` read-only, stress comparison N/A, and `*DLOAD` outside CLI support.