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.

No open issue requires formulation revision. NR-O01/NR-O02 are implementation-planning handoffs; NR-O03 belongs to 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

• Make NR-O03 orientation and unloaded-interior assumptions explicit without modifying the approved legacy 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.