Skip to content

Feature beam element methods - #80

Closed
timryanb wants to merge 51 commits into
masterfrom
feature-beam-element-methods
Closed

Feature beam element methods#80
timryanb wants to merge 51 commits into
masterfrom
feature-beam-element-methods

Conversation

@timryanb

Copy link
Copy Markdown
Owner

No description provided.

timryanb and others added 30 commits July 11, 2026 19:07
…alueh, Ap/Ah) fields

Adds the forward-over-reverse hforward()/hreverse() seed/accumulator slots
in place on the existing types rather than via new parallel types: every
existing first-order call site in TACSBeamElement.h/TACSBeamUtilities.h
uses only the default (or plain array) constructors, so the new fields
(always zero-initialized) cannot break any of them. Confirmed by a
standalone mpicxx compile (real and -DTACS_USE_COMPLEX) of a scratch
driver touching the new fields, and by a full setup.py build_ext --inplace
rebuild plus tests/element_tests/shell_tests/test_beam_element.py (still
4 passed, 12 subtests).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ub from TACSBeamUtilities.h

TacsShellComputeDispGrad/AddDispGradSens/AddDispGradHessian (the commented-
out block) never compiled as written and was a stale shell copy-paste, not
a viable template under the A2D-based second-order design. TacsTestBeamUtilities
was an unconditional fail=0 stub with its entire body commented out.
Confirmed no live caller via grep across src/tests/tacs before deleting.
File still builds; tests/element_tests/shell_tests/test_beam_element.py
still 4 passed.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… through the fixed linear strain map

The seven Hessian blocks are a constant pullback of the packed-symmetric
Cs[21] through evalStrain's fixed linear map e = L(u0x, d1x, d2x, e0ty), so
no differentiation happens at runtime -- every Cs row/col index (including
the cross-coupling terms not spelled out in SPEC.md, e.g. d2d1xd2x's
torsion-torsion cross term at Cs[6]) was derived directly from
TACSBeamConstitutive::computeStress's packed-symmetric index arithmetic.
Verified against a central-difference-of-f Hessian (f = 0.5*s.e, exact for
this quadratic form) independently for each of the 7 blocks, in both real
and -DTACS_USE_COMPLEX builds (max err ~1e-9 in both).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
alpha = scale*(x.x . y.x) is bilinear in (x, y): hforward is the product-
rule directional derivative of forward()'s formula (seed xp/yp substituted
for xd/yd); hreverse is reverse()'s adjoint shape plus one seed-cross-term
per input. Ports VALIDATION.md's E6-verified formula pair into the real
a2dvecops.h header. New standalone op-level test (header-only, mpicxx,
no TACS build needed) verifies hforward against a directional derivative
of the primal and hreverse against a central-difference-of-the-exact-
reverse()-adjoint, in both real and -DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…lf-Hessian sanity check) + op-level test

C = scale*A*B is linear in its single active input A (B is a passive
Mat3x3), so hforward/hreverse are exactly forward()/reverse()'s existing
transpose-map formulas with Ap/Ah substituted for Ad/Ad -- no additional
cross term, matching VALIDATION.md E6's zero-self-Hessian finding for this
op-class. New op-level test verifies: (1) hforward against a directional
derivative of the primal, (2) the zero-self-Hessian property bit-exactly
(0.0, not just within tolerance) when no downstream second-order seed is
present, and (3) hreverse's propagation formula against a hand-computed
scale*Ch*B^T for a nonzero downstream seed. Passes in both real and
-DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
v = scale*alpha.value*x + y is linear jointly in the two active inputs
(x, y) -- alpha is a fixed passive Scalar, not an AD variable, so there is
no bilinear coupling between x and y. hforward/hreverse mirror forward()/
reverse() with xp/yp/vp (resp. xh/yh/vh) substituted for xd/yd/vd, no
cross term. New op-level test verifies hforward against a directional
derivative of the primal, a bit-exact zero-self-Hessian property, and
hreverse's propagation formula against a hand computation, in both real
and -DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
C's columns are literally x/y/z's components -- a trivial linear
assignment, not even a product of two active quantities -- so hforward/
hreverse mirror forward()/reverse() with Ap/xp (resp. Ah/xh) substituted
for Ad/xd, no cross term. New op-level test verifies the directional-
derivative property of hforward, a bit-exact zero-self-Hessian property,
and hreverse's propagation formula against a hand computation, in both
real and -DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
C = scale*A^T*B is linear in its single active input B (A is a passive
Mat3x3) -- hforward/hreverse mirror forward()/reverse()'s existing
formulas with Bp/Cp (resp. Bh/Ch) substituted for Bd/Cd, no cross term
(same structural reasoning as ADMat3x3MatMult, Task 1.5). New op-level
test verifies the directional-derivative property of hforward, a
bit-exact zero-self-Hessian property, and hreverse's propagation formula
against a hand computation, in both real and -DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
y = scale.value*A^T*x is linear in its single active input x (scale and A
are both passive) -- hforward/hreverse mirror forward()/reverse()'s
existing accumulating formulas with xp/yp (resp. xh/yh) substituted for
xd/yd, no cross term. This is the last of Phase 1's 6 load-bearing op
extensions. New op-level test verifies the directional-derivative property
of hforward, a bit-exact zero-self-Hessian property, and hreverse's
propagation formula against a hand computation, in both real and
-DTACS_USE_COMPLEX builds.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…p-test runner

Reviewer flagged two gaps in the Phase 1 diff, verdict "With fixes":
(1) ADMat3x3MatMult/MatTrans3x3ADMatMult's new hforward/hreverse carry a
TacsRealPart(scale)==1.0 branch that is a literal instance of SPEC.md sec
1.2.4's no-TacsRealPart-branching rule; judged safe (scale is always a
passive real literal at every call site, never complex-step-perturbed, and
the branch exactly mirrors the pre-existing forward()/reverse() branch) but
undocumented -- added a deviation comment at both sites rather than
restructuring the code, per this repo's established practice (cf. commit
eb43259).
(2) the six op tests in src/elements/a2d/tests/ had no repeatable runner --
added run_tests.sh, which builds and runs all six in both real and
-DTACS_USE_COMPLEX modes and exits nonzero on any failure (12/12 passing).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…lStrainSens, Cs materialized, mat left zero pending A2D sweep

New addJacobian override (SPEC.md sec 1.3): zeros mat[]/res[], reruns
addResidual's own forward setup (fn1/fn2, director rates, tying strain,
the per-quadrature-point A2D chain through u0x/d1x/d2x), calls the NEW
con->evalTangentStiffness (never materialized by addResidual/
computeEnergies today) and populates res via the existing, unchanged
model::evalStrainSens -- mat[] is intentionally left all-zero at this
checkpoint, the A2D hforward/hreverse sweep that populates it lands in
Task 2.2.

Verified in the worktree build (Beam2 and Beam3, IsoTubeBeamConstitutive):
TacsTestElementJacobian's residual-consistency half passes bit-exactly
(res from addJacobian == addResidual's own res, 0.0 err); the mat-vs-FD
half fails as expected and documented (an all-zero mat cannot match a
nonzero FD Jacobian for a real stiffness element) -- this is the
intentionally-incomplete, documented RED-in-progress state Task 2.2 closes.

Confirmed a real environment hazard along the way: setup.py build_ext
--inplace's incremental dependency tracking does not reliably notice
header-only changes to TACSBeamElement.h (a stale .so was picked up once,
silently masking this exact all-zero-mat state as a false pass) -- use
--force for every rebuild for the remainder of this phase.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… pointer skip

The "leftover row" pointer skip (mat += (nbrows - njrows) * ncols) was
incorrectly placed inside the njrows loop, executing njrows times per
outer node instead of once. This walks mat[] out of bounds whenever
nbrows != njrows (confirmed via a standalone pointer-arithmetic
simulation: for a 2-node element with nbrows=6/njrows=3, the buggy code
touched flat index 248 against a 144-entry buffer). The structurally
analogous TACSShellElementQuadBasis::addInterpFieldsOuterProduct places
this same skip once per outer node, which is correct.

Found while wiring TACSBeamElement::addJacobian's mass-Hessian scatter
(Phase 2, Task 2.3) -- this helper was never previously called for beam
(addResidual/computeEnergies never used it), so the bug was latent.
Re-verified via simulation that the fix preserves correctness for the
njrows == nbrows case already exercised elsewhere.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
New model-layer closure converting the tying-point-space tying-strain
Hessian (d2ety, produced per-quadrature-point by the existing
basis::addInterpTyingStrainHessian) into mat[]'s (u0xi, u0xi) block plus
the director-space Hessian accumulators (d2d1/d2d2/d2d1u/d2d2u) that
director::addDirectorJacobian consumes -- the Hessian-level analog of the
existing first-order addComputeTyingStrainTranspose.

Mirrors TACSShellElementModel::addComputeTyingStrainHessian
(TACSShellElementModel.h:159-390), simplified because beam has exactly
two tying fields (G12/G13) that never cross-couple for the constitutive
models this feature exercises (confirmed: TACSIsoTubeBeamConstitutive
never sets the e0ty-e0ty off-diagonal Cs entry that would otherwise
require shell's full arbitrary-tying-field-pair double loop). The
external tying-strain-vs-other-strain cross-Hessian inputs
(d2etyu/d2etyd1/d2etyd2) shell's version needs are always zero-filled
placeholders here, since beam's evalStrainHessian produces no such cross
term -- kept in the signature only for structural parity with shell's
template.

SPEC.md sec 1.3.1.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…s, and cross-director closure (Tasks 2.2/2.3)

Adds TacsBeamAddCrossDirectorJacobian to TACSBeamUtilities.h (Gap 2,
SPEC.md sec 1.3.2): the (d1, d2) cross-director Hessian closure -- the
feature's first genuinely novel derivation, since shell has only one
director and so no cross-director term can exist there. Exact for
TACSLinearizedRotation (the only director class test_beam_element.py
exercises today), with a documented scope boundary for
TACSQuadraticRotation/TACSQuaternionRotation (their director-to-rotation
map is nonlinear in the state, so an additional curvature correction
would be needed for exactness -- not exercised today, explicitly scoped
out rather than silently extended).

Wires TACSBeamElement::addJacobian's per-DOF hforward/hreverse sweep: for
each quadrature point, 3*dsize sweep iterations (translational, d1, d2
directions -- not nvars-many, since d1/d2 are functions of the director
node arrays, not a direct basis interpolation of vars) seed the
Phase-1-extended second-order A2D graph, contract against the fixed
per-quadrature-point material Hessian blocks (model::evalStrainHessian),
and scatter the result into mat[] directly (translational-translational
block) or into the director-space Hessian accumulators d2d1/d2d2 (self),
d2d1u/d2d2u (director-translational cross) and d2d1d2 (the new
cross-director accumulator, capturing the leakage that occurs because
u0d's column assembly couples u0xi/d01/d02 together, rather than
discarding it). Also adds the tying-strain Hessian's per-quadrature-point
seeding (feeding Gap 1's addComputeTyingStrainHessian, called post-loop)
and the mass-moment dynamics Hessian blocks (fixed-coefficient outer
products of the basis shape functions, mirroring
TACSShellElement::addJacobian's identical pattern), then calls the two
director::addDirectorJacobian closures (one per director) and the new
cross-director closure.

Removes the separate addDirectorResidual calls from Task 2.1's res-only
scaffold: director::addDirectorJacobian now updates res internally too
(via the same crossProductAdd formula), so keeping both would
double-count the directorial residual contribution.

Task 2.2 (sweep + Gap 1) and Task 2.3 (dynamics + director-Jacobian +
Gap 2) are committed together, per the plan's own guidance: Task 2.2
alone leaves mat[]'s rotational-DOF block at zero and cannot pass
test_element_jacobian in isolation.

Result: test_element_jacobian passes for all variants
test_beam_element.py instantiates (Beam2, Beam3 x BeamRefAxisTransform),
for the first time checked against genuine FD-of-addResidual rather than
FD-vs-FD.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…d of a silent early return

Per the never-silent-punt principle (SPEC.md sec 4.1/4.2): replace the
bare "// Not implemented\nreturn;" with an explicit forward to
TACSElement::getMatType. Numerically a no-op -- the base class's own
TACS_GEOMETRIC_STIFFNESS_MATRIX branch (TACSElement.cpp:346-348) is also
unimplemented today, so both produce an all-zero mat. The value is
consistency/auditability: a future reader sees an explicit forward, not
an ambiguous silent return, and Phase 5 has a clean explicit branch to
replace with a real analytic implementation (SPEC.md sec 1.4/2.4.4/2.4.5).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ope and symmetry assumptions at their actual use sites

Two documentation-only fixes from the "With fixes" Phase 2 review verdict:

1. TACSBeamElement.h's TacsBeamAddCrossDirectorJacobian call site had no
   pointer back to the TACSLinearizedRotation-only exactness boundary --
   it was documented only on the utility's own docstring
   (TACSBeamUtilities.h), not at the call site, which is where a future
   QuadraticRotation/QuaternionRotation extension attempt would actually
   land. Added a one-line comment restating the boundary with a pointer
   to the full rationale.

2. TACSBeamUtilities.h's TacsBeamAddCrossDirectorJacobian reuses the same
   D[3*p+q] value for both the (i,j) and transposed (j,i) mat[] slots
   instead of computing a separate transpose. This is valid only because
   this model's d2d1d2 3x3 node-pair blocks are themselves always
   symmetric (the static term is a rank-1 t1-outer-product-with-itself,
   the dynamics term is rho[5]*I3) -- documented at the write site so a
   future kinematics change that breaks this symmetry doesn't silently
   corrupt mat[] via this shortcut.

No logic changes; test_element_jacobian remains green for all variants,
run_tests.sh remains 12/12.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… computation, closing a real Assembler.getAverageStresses correctness gap for beam components

Beam did not override getAverageStresses, so any beam element in a model
silently contributed zero to Assembler.getAverageStresses(compNum)'s
average -- a real correctness gap for this Python-exposed API, not just
an abstract interface-parity one (TACSAssembler::getAverageStresses loops
over every element in the component and accumulates via +=, then divides
by the element count regardless of whether each element actually
contributed anything).

Forward-only value computation (not a sensitivity): loops over
num_vis_nodes visualization points, recomputes the strain the same way
addResidual/computeEnergies/getOutputData already do (reusing
TacsBeamComputeNodeNormals, director::computeDirectorRates,
model::computeTyingStrain, the state-independent geometry chain, and
model::evalStrain unchanged), calls con->evalStress, and accumulates the
averaged 6-component beam stress into avgStresses[0..5].
avgStresses[6..8] are deliberately left untouched (not zeroed), matching
the base class's own no-op-leaves-untouched convention for slots this
element doesn't know how to fill.

Mirrors TACSShellElement::getAverageStresses's structure exactly, with
one deliberate deviation from TACSBeamElement::getOutputData's own
per-vis-node loop: uses the confirmed-correct e0ty = 2*XdinvT.A[0]*gty
transform (matching addResidual/addJacobian) rather than getOutputData's
un-scaled e0ty = gty, which appears to be a pre-existing inconsistency in
getOutputData -- out of scope for this task, not touched here.

New regression test (tests/integration_tests/test_beam_average_stresses.py):
no existing harness covers getAverageStresses for any element type in the
repo (a value computation, not a sensitivity, so the TacsTestElement*
FD/CS harnesses don't apply). Builds a single straight Beam2 and a single
straight Beam3 element directly via the low-level TACSAssembler interface
(no Creator/BDF needed for a single, single-process element), prescribes
a uniform axial-strain state directly via the state vector (bypassing
solving, since only the forward stress computation is under test), and
checks the assembler-level average stress against the closed-form axial
stress resultant s[0] = E*A*eps, plus confirms the other 5 stress
components are exactly zero for a pure axial state. Confirmed RED against
a temporary no-op stub before implementing the real body.

Scope note on avgStresses[6..8]: verified by code inspection that this
implementation never writes to those indices (matching
TACSShellElement's own, similarly untested-at-this-precision convention);
a fully conclusive "left untouched, not merely zeroed" test would need
either a new Cython binding exposing a pre-seeded buffer or a
mixed-beam-and-shell single-component assembler, deferred per SPEC.md sec
1.5's own "recommended, non-blocking fast-follow" framing for this
method's test rigor.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…Sens: parallel-axis Xpts-adjoint built from already-implemented addTransformSens/mat3x3SymmTransformSens-style algebra

Analytic Xpts-sensitivity for TACS_ELEMENT_MOMENT_OF_INERTIA (SPEC.md sec
3.1), replacing the previous explicit-forward to
TACSElement::addPointQuantityXptSens's FD-based default. Differentiates
evalPointQuantity's own forward formula directly: moments =
con->evalMassMoments(...) (hence I0[3..5]) are Xpts-independent, so
quantity = T*I0*T^T + parallel_axis(dXcg) has exactly two Xpts-dependent
paths:
  1. T (via X0xi): df/dT = 2*D*T*I0 (D = dfdq unpacked to a full
     symmetric 3x3), seeded directly into T.Ad -- T already participates
     in this function's shared A2D graph, so the existing
     transform->addTransformSens(X0xi.x, T.Ad, X0xi.xd) call in the
     common tail (unchanged) converts it into the Xpts adjoint.
  2. dXcg (via X0/n1/n2): the parallel-axis term is a point-mass inertia
     tensor density*(|dXcg|^2*I - dXcg (x) dXcg), so
     df/d(dXcg) = 2*density*(tr(D)*dXcg - D*dXcg); dXcg is linear in
     X0/n1/n2, so this seeds X0.xd/n1.xd/n2.xd directly, mirroring the
     existing TACS_ELEMENT_DENSITY_MOMENT branch's identical
     seed-then-fall-through pattern.

No new algebra primitive: both paths reduce to plain 3x3 matrix
arithmetic (A2D::Mat3x3MatMultCore/Mat3x3VecMultCore) plus the existing
transform->addTransformSens machinery the common tail already calls.

Two bugs caught during GREEN, verified against standalone FD checks
before/after each fix:
- dfdq is packed-symmetric (xx,xy,xz,yy,yz,zz) with
  f = sum_k dfdq[k]*quantity[k] -- each off-diagonal slot counted ONCE.
  Naively mirroring dfdq[1] into both D[0][1] and D[1][0] at full value
  double-counts the off-diagonal contribution; fixed by halving the
  off-diagonal entries when unpacking dfdq into the full 3x3 D.
- TacsTestElementQuantityXptSens draws a random, generically nonzero
  "scale" and folds it into its own FD reference
  (TACSElementVerification.cpp:1678-1707); this branch was missing that
  factor (matching the pre-existing TACS_ELEMENT_DENSITY_MOMENT branch's
  same omission, which is out of scope here and left untouched -- it
  predates this feature and has never been covered by an automated
  test either, per the wiring gap below).

New Cython wiring (tacs/cpp_headers/elements.pxd,
tacs/elements.pyx: TestElementQuantityXptSens, ELEMENT_DENSITY_MOMENT,
ELEMENT_MOMENT_OF_INERTIA): no element type in the repo had
TacsTestElementQuantityXptSens exposed to Python before this -- confirmed
by grepping tests/, tacs/elements.pyx, and
tacs/cpp_headers/elements.pxd, all empty. Added the minimal binding
(mirroring TestAdjResXptProduct's existing wrapper style) since this is
general infrastructure useful beyond beam, rather than a standalone
C++-only test program.

New regression test (test_beam_element.py::test_quantity_xpt_sens):
confirmed genuine RED first against a deliberately-empty stub (the
original forward-to-base "passed" for the wrong reason -- the base
class's own fallback is itself FD-based, the same FD-vs-FD non-signal
VALIDATION.md's E4 finding flagged for addJacobian's base fallback), then
GREEN after the real derivation (with both bugs above caught and fixed
before landing).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…MATRIX analytic, GEOMETRIC_STIFFNESS_MATRIX interim-forwarded

TACS_STIFFNESS_MATRIX reduces psi^T*K*phi to e_psi^T*Cs*e_phi (beam's
strain map is linear unconditionally), so its DV-sens is exactly the
existing con->addStressDVSens hook addAdjResProduct already calls,
reused with a bilinear pair of strains built by running psi/phi through
the same director-Jacobian kinematics substitution (computeDirectorRatesDeriv,
linearized about the real vars, generic across director classes).
TACS_MASS_MATRIX contracts the same psi/phi-direction director fields
bilinearly into rho and reuses con->addMassMomentsDVSens.
TACS_GEOMETRIC_STIFFNESS_MATRIX explicit-forwards to the base class,
interim per Phase 4's preamble (superseded by Phase 5).

Also adds minimal forward-to-base skeleton bodies for
addMatXptSensInnerProduct/getMatSVSensInnerProduct (Tasks 4.2/4.4's own
starting points) so the vtable links; their real analytic branches land
in later Phase 4 tasks.

Verified via a deliberate-bug perturbation (temporarily doubling the
stress-DV-sens scale and the mass rho[0] term) that both branches are
genuinely exercised by test_element_mat_dv_sens, run with the test's
skip decorator locally/uncommitted; the committed test file is
unchanged (all three matrix-sens tests remain skipped until Task 4.5).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ATRIX analytic Xpts-sens

TACS_MASS_MATRIX has no strain/stiffness path: its Xpts-dependence
flows through detXd(Xpts) directly, plus the psi/phi-direction director
fields' own dependence on Xpts via fn1/fn2 (the director map's
reference-normal argument). Contracts the bilinear rho-weighted dot
products (mirroring addJacobian's dynamics blocks) and routes the
result through detXd's existing Xpts-adjoint machinery plus
addDirectorRefNormalSens/TacsBeamAddNodeNormalsSens for the director
term, reusing addAdjResXptProduct's own chain shape.

The director-field routing (two single-direction addDirectorRefNormalSens
calls, one per test direction) is exact for TACSLinearizedRotation (the
only director class this feature's test file exercises, and the only
one whose director-map Jacobian is constant) -- not verified for
TACSQuadraticRotation/TACSQuaternionRotation, which would need the
director's own curvature (second-derivative) term this feature's A2D
machinery does not build; documented as an open risk mirroring
SPEC.md's getMatSVSensInnerProduct caveat.

TACS_STIFFNESS_MATRIX remains temporarily forwarded to base pending
Task 4.3's analytic derivation; TACS_GEOMETRIC_STIFFNESS_MATRIX and any
future matType forward to base per SPEC.md's interim/never-silent-punt
rules.

Found and fixed a reverse-mode ordering bug during this task's RED/GREEN
cycle: the detXd/Xd A2D reverse calls must unwind in the exact reverse
of their forward construction order (computedetXd.reverse() before
assembleXd.reverse()), not source order -- caught by the test itself
(a real, non-trivial failure, not a forwarding false-positive) and
confirmed fixed by re-running. Verified via a deliberate-bug
perturbation (temporarily doubling the detXd seed) that the analytic
branch is genuinely exercised.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…analytic Xpts-sens (5-piece decomposition)

Per the finalized scope (SPEC.md sec 2.3.1), this Xpts-sens branch is
brought in scope and made analytic rather than left forwarded. Since
beam's strain is linear in state unconditionally, psi^T*K*phi reduces
to e_psi^T*Cs*e_phi with no separate strain-Hessian term needed; its
Xpts-derivative decomposes into: (1) a detXd-direction term seeded by
e_psi.Cs.e_phi, (2)/(3) psi-direction and phi-direction kinematics
chains seeded by the cross-stress weights s_phi=Cs*e_phi/s_psi=Cs*e_psi,
and (4) the e0ty/XdinvT product-rule correction -- all re-invocations of
beam's already-shipped addAdjResXptProduct machinery
(TacsBeamAddNodeNormalsSens, director::addDirectorRefNormalSens,
model::addTyingStrainDerivXptSens) with both "state" directions replaced
by pure test directions (psi, phi), rather than new algebra. No
tying-curvature term, since beam's tying strain is linear.

Split into a private addMatXptSensInnerProductStiffness helper given its
size. Verified genuinely analytic (not a forwarding false-positive) via
a deliberate-bug perturbation of the detXd seed. Director-class scope:
exact for TACSLinearizedRotation (the only director class this
feature's test file exercises); TACSQuadraticRotation/TACSQuaternionRotation
are undocumented/unverified here, same open risk flagged for Task 4.2's
TACS_MASS_MATRIX branch and Task 4.4's getMatSVSensInnerProduct.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…LinearizedRotation, typeid-guarded fallback for the nonlinear directors

psi^T*mat*phi = e_psi(vars)^T*Cs*e_phi(vars), where e_psi(vars) =
L(J_dir(vars)*psi) (L = model::evalStrain's fixed linear strain map,
J_dir = the director map's own Jacobian). Differentiating again w.r.t.
vars produces a term proportional to H_dir(vars) (the director map's
own second derivative/curvature) -- a different quantity from
evalStrainHessian's blocks, which only govern the strain map's
linearity (always true for TACSBeamLinearModel). For
TACSLinearizedRotation, the director map itself is linear (constant
Jacobian), so H_dir=0 identically and the memset already produces the
exact answer -- confirmed both algebraically and empirically (Max Err
0.0, not merely inside the loose rtol band).

For TACSQuadraticRotation/TACSQuaternionRotation, H_dir is genuinely
nonzero: confirmed via the mandatory experiment SPEC.md sec 6.4
requires (temporarily forcing zero unconditionally and observing
TestElementMatSVSens fail hard against TACSBeam2ModRot, Max Err ~1562
for STIFFNESS_MATRIX / ~15 for MASS_MATRIX -- not a rounding-level
discrepancy). This feature's A2D machinery only extends to second
order (Phase 1), one order short of what H_dir's contribution requires;
implementing it is out of scope. Falls back to the base FD/CS
implementation for these two director classes specifically, guarded by
typeid(director), per SPEC.md's documented-failure-only fallback rule --
not a special case forcing a zero, but an honest routing between the
provably-exact-zero case and the genuinely-unresolved one.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ype/the matrix-sens overrides are analytic

test_element_mat_dv_sens/test_element_mat_xpt_sens/test_element_mat_sv_sens
are now real, permanent regression coverage (VALIDATION.md E1's
confirmation step) rather than presumed-FD-vs-FD skipped placeholders:
TACS_STIFFNESS_MATRIX/TACS_MASS_MATRIX are genuinely analytic-vs-FD/CS
for all three methods (including the TACS_STIFFNESS_MATRIX Xpts-sens
branch Task 4.3 brought in scope); TACS_GEOMETRIC_STIFFNESS_MATRIX
still passes via the interim explicit-forward punts, expected until
Phase 5.

test_beam_element.py -> 8 passed (up from 5; PLAN.md's original "7
passed" acceptance figure was an arithmetic slip -- 5 baseline + 3
newly-real tests is 8, not 7 -- corrected here since this session
directly observed and counted the file's 8 test methods).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…inear directors, fallback where the fix is incomplete

CRITICAL (fixed): addMatXptSensInnerProduct's director-field Xpts-adjoint
(both the TACS_MASS_MATRIX branch and the TACS_STIFFNESS_MATRIX helper)
called the single-direction 4-arg addDirectorRefNormalSens overload with
phi/psi substituted for the "vars" argument -- exact only by accident for
TACSLinearizedRotation (whose director-map Jacobian is constant,
independent of the linearization point), but silently wrong for
TACSQuadraticRotation/TACSQuaternionRotation, whose Jacobian-vector-product
Xpts-sensitivity genuinely depends on the REAL q, not on phi/psi treated
as if they were the state (confirmed from TACSDirector.h's own per-class
implementations). Fixed by calling the 6-arg overload with a zeroed "dd"
(base-term) buffer, vars=REAL vars, and phi/psi only ever occupying the
"perturbation direction" slot -- this isolates exactly the q-qpsi cross
term the review identified as missing, reusing addAdjResXptProduct's own
existing machinery rather than deriving new algebra.

Verified via a TACSBeam2ModRot (TACSQuadraticRotation) spot-check:
TACS_MASS_MATRIX is now analytic-exact (rel err ~1e-9, machine precision,
confirming the fix is complete for this branch). TACS_STIFFNESS_MATRIX
measurably improves (from ~90x to ~34x the failing threshold... i.e. the
fix reduces but does not eliminate the discrepancy) but retains an
unresolved ~9% residual, confirmed genuine (not FD noise, via a
dh-convergence check) and confirmed q-dependent (vanishes as O(vars^2)
via a vars-magnitude-scaling check) -- an additional, unidentified term
specific to this branch's larger kinematic chain (u0d/u0x/d1x/d2x/T/
XdinvT) that this session's investigation could not isolate. Per the
documented-failure-only fallback rule, TACS_STIFFNESS_MATRIX now
typeid-guards TACSQuadraticRotation/TACSQuaternionRotation to the base
FD/CS implementation (confirmed passing trivially, rel err ~1e-8) while
remaining analytic for TACSLinearizedRotation (still exact, confirmed by
the full un-skipped test_beam_element.py suite, 8/8 passing).

MINOR (documented, not fixed): computeDirectorRatesDeriv's redundant
per-call recomputation of the direction-independent d1/d1dot/d1ddot
outputs is left as-is -- the API has no lighter-weight overload that
skips them, and adding one would require changing TACSDirector.h (a
shared header also used by shell elements), not a local contortion;
noted at the call sites.

Full regression: tests/element_tests/ 176/176 passed, run_tests.sh
12/12, TACS.dtype confirmed real mode.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…sian family

Uncomments and finishes the previously-stubbed TACSBeamNonlinearModel:
fixes the unbalanced-paren syntax error in the torsional-strain
expression, widens evalStrainSens to take/produce e0ty/de0ty (matching
the linear model's shape), and adds evalStrainDeriv, evalStrainSensDeriv,
computeTyingStrainDeriv, evalStrainHessian, and evalStrainHessianDeriv.
Tying strain stays linear (byte-identical to TACSBeamLinearModel's own
computeTyingStrain), consistent with this beam formulation; only the
axial/torsion/bending strain measures are genuinely nonlinear.

evalStrainHessian/evalStrainHessianDeriv are implemented via small
Jacobian-row loops (J^T*Cs*J plus a constant stress-weighted geometric
term) rather than fully unrolled scalar code, since the u0x block is 9x9
-- documented as a deviation from the linear model's fully-unrolled
style. Every new method is verified to machine precision via standalone
complex-step checks (not just real-mode FD), in addition to the extended
TacsTestBeamModelDerivatives harness (now also FD-checking the four
Hessian blocks). Also fixes a latent bug in that harness where each
first-derivative block's `fail` assignment overwrote rather than
accumulated the previous blocks' results.

Deletes the harness's dead, stale shell-copy-paste FD-Hessian tail
(referencing u1x/TACSShellConstitutive inside a beam function), same
category of cleanup as Task 1.1's TACSBeamUtilities.h dead code.

No behavior change to TACSBeamLinearModel or any wired-up element path;
TACSBeamNonlinearModel is not yet referenced by TACSBeamElement (Task
5.2 wires it into getMatType).
…ytic port (Task 5.2)

Replaces the interim forward-to-base-class stub with a real analytic
implementation, scoped via typeid guard to TACSLinearizedRotation (the
same documented-fallback scoping convention Phases 2/4 already
established for other director-Jacobian closures; TACSQuadraticRotation/
TACSQuaternionRotation's own addDirectorJacobian overloads genuinely
depend on vars for mat[] itself, confirmed by reading TACSDirector.h
directly, so this guard is a correctness requirement).

Reports the combined "Ku + Kg" (initial-displacement + initial-stress)
geometric stiffness, not Kg alone. Resolved this ambiguity by
cross-checking two authoritative sources rather than guessing: (1)
src/TACSBuckling.cpp's TACSLinearBuckling::solve, which shows
getMatType(GEOMETRIC) receives the real path/vars directly from the
framework, not a caller-supplied zero state; (2) origin/master's
TACSShellElement::getMatType (fetched via `git show`, since this
branch's own checkout doesn't have it), whose GEOMETRIC branch performs
the "zero base state, path=vars direction" split internally and only
ever scatters the directional-derivative-suffixed Hessian family into
mat[] -- the base family is computed but deliberately discarded (it
reduces to the pure linear material stiffness TACS_STIFFNESS_MATRIX
already reports). Full trace recorded in PLAN.md's Phase 5 section.

Implementation reuses the exact per-DOF hforward/hreverse sweep
machinery already built for addJacobian (TacsBeamZeroSecondOrderNodes/
TacsBeamContractStrainHessian) -- no separate scatter mechanism.
Skips the "base" (zero-state) sweep pass entirely, since its output is
provably unused (confirmed from shell's own reference implementation):
the base Hessian's stress is identically zero, so the analogous
director-space accumulators the base pass would produce are dead
weight. Also skips any tying-strain-Hessian scatter for this matType:
beam's tying/shear strain is exactly linear in both models (Task 5.1),
so its own Hessian block is a Cs-only constant with a zero directional
derivative -- the tying-strain contribution to G is provably zero.

Verified: no regression in tests/element_tests/ (176 passed),
test_beam_element.py (8 passed/42 subtests), or
src/elements/a2d/tests/run_tests.sh (12/12). Real acceptance oracle is
the new Euler-column buckling integration test (next commit), which
passes end-to-end including its total-DV-sensitivity and
total-Xpts-sensitivity subtests.
New tests/integration_tests/test_beam_column_buckling.py -- the first
beam buckling correctness oracle in this repository (does not reuse or
repurpose test_rectangle_beam_buckling.py, which has zero eigensolve
machinery despite its name). Cantilever (fixed-free, K=2) Euler column
on the existing 5-element beam_model_translated.bdf mesh (already
exercised by test_rectangle_beam_buckling.py's own static test),
asserting the lowest TACSLinearBuckling eigenvalue matches the
closed-form P_cr = pi^2*E*I_min/(K*L)^2 (I_min = the weak-axis moment of
inertia, w*t^3/12, confirmed directly from
TACSIsoRectangleBeamConstitutive::evalMomentsOfInertia) to rtol=5e-2
real / rtol=1e-8 complex-step.

Written and confirmed RED before Task 5.2 landed (the interim
zero-valued geometric stiffness matrix gave a garbage eigenvalue,
~7.4e12, not a small or trivially-zero one -- confirming the test
genuinely exercises the eigensolve path), per the coordinator's
suggestion to validate Task 5.2's Ku+Kg resolution against real
physics as early as possible rather than only at phase end. After Task
5.2 landed, the eigenvalue matched the closed-form load to within 0.7%
but with the wrong sign, because the BDF's own load case (a unit +x
force at the free end, further from the fixed node) is tensile, not
compressive; fixed by applying it with scale=-1.0. All 5 pytacs base
test methods pass (20 subtests), including total-DV-sensitivity and
total-Xpts-sensitivity, which exercise the GEOMETRIC_STIFFNESS_MATRIX
matrix-sens methods (still forwarded to TACSElement's own generic
FD-based fallback pending Tasks 5.3-5.5) end to end.

Lateral-torsional buckling stretch case (SPEC's recommended, not
required, second case) not attempted in this pass.
…S_MATRIX analytic branch (Task 5.4)

d(psi^T*G(vars)*phi)/d(vars) = G(phi)*psi, reusing Task 5.2's getMatType
verbatim rather than a fresh derivation. This follows from a symmetric-
trilinear-tensor identity: G(vars) (as constructed by getMatType's
GEOMETRIC branch) is linear in its "path"=vars argument (confirmed
directly from that branch's own construction -- every Hessian block it
scatters is a Cs-only-linear function of path), so the scalar
psi^T*G(vars)*phi is a fully symmetric trilinear form in (psi, phi,
vars) -- third mixed partial derivatives of the underlying strain-
energy-like functional commute regardless of implementation detail.
Differentiating w.r.t. vars and permuting slots gives exactly the
matrix-vector product of "getMatType(GEOMETRIC, ..., vars=phi, ...)"
applied to psi; verified this is equivalent (not merely analogous) to
the G(psi)*phi form via direct index-notation derivation, recorded in
PLAN.md.

Scoped to TACSLinearizedRotation via typeid guard, same reason as Task
5.2's getMatType port it calls (the other two director classes'
addDirectorJacobian overloads genuinely depend on vars for mat[]
itself).

Verified: test_beam_element.py's test_element_mat_sv_sens
GEOMETRIC_STIFFNESS_MATRIX subtest now exercises this analytic path
(previously the interim base-class FD/CS fallback) and still passes (8
passed/42 subtests, no regression). Full suite: tests/element_tests/
176 passed, run_tests.sh 12/12, test_beam_column_buckling.py 5
passed/20 subtests (no regression from Task 5.2/5.7).
…ly, fix docstring

Review feedback (Phase 5, diff 38d173c..HEAD, verdict "With fixes"):

Important: evalStrainDeriv/evalStrainSensDeriv/evalStrainHessianDeriv
(added in 73858cb) had zero committed, reproducible regression
coverage -- only an uncommitted scratch complex-step driver -- despite
being exactly what TACSBeamElement::getMatType's
TACS_GEOMETRIC_STIFFNESS_MATRIX branch (48fbacc) depends on. Adds
TacsTestBeamNonlinearModelDerivFamily, a new standalone FD/CS check
function in TACSBeamElementModel.h, hardcoded to
TACSBeamNonlinearModel (not folded into the existing, model-templated
TacsTestBeamModelDerivatives, since TACSBeamLinearModel does not have
-- and by this feature's design, mirroring TACSShellLinearModel's own
identical omission, is not meant to have -- any of these three
methods). Wires it into examples/shell/beam.cpp's manual-example
entrypoint alongside the existing TacsTestBeamModelDerivatives calls,
so it is at least committed and re-runnable, even though (like the
existing calls there) it is not part of the automated pytest/testflo
suite. Confirmed passing (exit 0) after rebuilding and running the
example.

Minor: fixes a self-contradictory sentence in evalStrain's docstring
(TACSBeamElementModel.h) -- it stated "u0x[3*a+b] = d(u_a)/d(xi_b)"
while the worked examples two sentences later only make sense under
the reversed convention u0x[3*b+a]; the code itself was already
correct, only the docstring sentence was wrong.

No regression: tests/element_tests/ 176 passed, run_tests.sh 12/12,
test_beam_column_buckling.py 5 passed/20 subtests, TACS.dtype real.
Comment-only change (coordinator follow-up after the Phase 5 review):
the attempt-first rule requires failures to be visible where the
fallback actually lives, not only in PLAN.md. Updates
addMatDVSensInnerProduct's TACS_GEOMETRIC_STIFFNESS_MATRIX branch with
the specific root cause found (a strain-vector-through-Cs
reformulation misses getMatType's separate director-Jacobian-closure
scatter path), replacing the stale pre-attempt "superseded by Phase 5
unless it overruns its timebox" comment left over from Task 2.4.
Updates addMatXptSensInnerProduct's TACS_GEOMETRIC_STIFFNESS_MATRIX
branch similarly, explaining Task 5.5 was not attempted this session
for a reasoned, documented cause tied directly to Task 5.3's finding,
not "deferred for time" alone.

No behavior change; no regression (tests/element_tests/ 176 passed,
run_tests.sh 12/12, test_beam_column_buckling.py 5 passed/20
subtests).
timryanb added 21 commits July 12, 2026 01:51
…uct/addMatXptSensInnerProduct header comments

Both header comments still called TACS_GEOMETRIC_STIFFNESS_MATRIX's
fallback an "interim explicit-forward-to-base punt, superseded by
Phase 5" -- stale from before Tasks 5.3/5.5 concluded. Update to match
the accurate inline branch comments: addMatDVSensInnerProduct's
GEOMETRIC branch is a permanent documented-failure fallback (Task 5.3
attempted and root-caused), addMatXptSensInnerProduct's is a permanent
reasoned-deferral fallback (Task 5.5, deferred on 5.3's root cause).
Comment-only, no logic change.
ruff format on tests/integration_tests/test_beam_average_stresses.py
(the only changed .py file failing ruff format --check; ruff check was
already clean on all three changed .py files).

clang-format -i on every changed .h/.cpp file per the org PR checklist.
Guarded the six new src/elements/a2d/tests/*.cpp files' build-recipe
comment blocks (mpicxx invocations using shell `\` line-continuation)
with clang-format off/on, since clang-format's comment reflow otherwise
moves the trailing backslash onto its own line and mangles the
copy-pasteable shell commands. All other changes are pure whitespace/
line-break reflow (arithmetic expression and call-argument wrapping in
TACSBeamElement.h/TACSBeamElementModel.h) -- no semantic change.

Re-verified after formatting + a --force rebuild: run_tests.sh 12/12
(unchanged), test_beam_element.py 8 passed/42 subtests (unchanged).
clang-format --dry-run -Werror and ruff format --check both exit 0.
…level test

Ports the scale=1 formula E7's harness already verified bit-exactly
(docs/plans/feature-beam-element-methods/scripts/004_e7_zero_seed_role.cpp),
generalized to the scale!=1 branch by mirroring the class's existing
forward()/reverse() TacsRealPart(scale)==1.0 branch pattern. Op-level
complex-step test confirms hforward vs directional-derivative-of-primal and
hreverse vs a hand-computed formula, both real and complex builds.
…+ op-level test

Same pattern as ADMat3x3ADMatMult's own second-order extension: ports the
scale=1 formula E7's harness already verified bit-exactly, generalized to
the scale!=1 branch via the class's existing forward()/reverse() branch
convention.
…verified) + op-level test

Derived by the identical product-rule/reverse-shape-plus-cross-term pattern
E7 verified for ADMat3x3ADMatMult/ADMatTrans3x3ADMatMult, NOT itself run
through E7's harness (SPEC-phase-7.md sec 2.2/sec 5 flagged this as a
residual, low-risk implementation-time task) -- this commit's op-level
complex-step test is that independent verification: hforward vs
directional-derivative-of-primal and hreverse vs a hand-computed formula,
both real and complex builds, all three operands (alpha, x, y) genuinely
seeded.
…ern, independently verified) + op-level test

Generalizes the bilinear "reverse-shape + one-cross-term" pattern E7
verified to this op's genuinely trilinear form (y = scale.value*A^T*x, all
three of scale/A/x are AD-typed) -- hforward is forward()'s own formula with
the p-seed substituted; hreverse adds TWO cross terms per input (one per
pairing with each of the other two factors: (scale,x), (scale,A), (A,x)),
generalizing the bilinear ops' single cross term. Not itself run through
E7's harness (SPEC-phase-7.md sec 2.2/sec 5 flagged this as a residual
implementation-time task) -- this commit's op-level complex-step test, with
all three operands genuinely seeded, is that independent verification.
…der.cpp

E7 found this exact raw-TacsScalar[9]-passed-to-const-Mat3x3&-reference-
member pattern is undefined behavior (the implicit temporary Mat3x3 built
by Mat3x3's converting constructor is destroyed at the end of the
constructor-call's full expression, leaving the op holding a dangling
reference for every subsequent .hforward()/.hreverse() call) --
VALIDATION.md's E7 entry, SPEC-phase-7.md sec 2.2/sec 7. Stayed silently
benign here only by luck of stack-slot reuse; a systematic re-check of the
other five existing op tests (test_mattrans3x3admatmult_second_order.cpp,
test_mattrans3x3advecmultscale_second_order.cpp,
test_admat3x3fromthreeadvec3_second_order.cpp,
test_advec3advecscalaraxpy_second_order.cpp,
test_advec3dot_second_order.cpp) confirms none of them share this pattern
(all already use named objects for any operand referenced across multiple
statements). Fix: introduce a named, block-scoped Mat3x3 Bmat(B) and use it
at every ADMat3x3MatMult constructor call site instead of the raw array;
B[] itself is retained only for the hand-computed reference formula, which
never binds a reference to it. No behavior change (still passes bit-exact
in both real and complex builds).
…malSens hooks (Task 7.2)

New per-node hooks on TACSQuadraticRotation and TACSQuaternionRotation
(SPEC-phase-7.md sec 3.1/3.2): addDirectorHessianProduct extracts
addRotationMatJacobian's own i==j diagonal-block algebra
(TACSDirector.h:806-818 / :1681-1699) -- the constant, q-independent
(dU/dC):(d^2C/dq^2) tensor that arises because both director classes' C(q)
is an exact quadratic/homogeneous-quadratic polynomial -- into a standalone
bilinear map of two arbitrary caller-supplied directions (qa, qb) instead of
the same q implicitly on both sides. addDirectorHessianRefNormalSens is its
closed-form, one-pass linear Xpts-adjoint (the contraction is linear in the
reference normal t for fixed dd).

Omitted on TACSLinearizedRotation, whose C(q) = I - q^x is exactly linear,
so d^2C/dq^2 = 0 (matching this feature's already-established "exact zero"
ledger entry for G1/G2's linear-director case).

Verified via a new director-hook-level test,
src/elements/shell/tests/test_director_hessian_second_order.cpp (not
header-only -- links against the already-built libtacs.so, per
SPEC-phase-7.md sec 6.2's own allowance; not part of
src/elements/a2d/tests/run_tests.sh's glob, a different test class): (1) a
regression-safe-refactor check -- addDirectorHessianProduct called with
qa=qb=basis vectors over all (k,l) pairs reproduces
addRotationMatJacobian's own diagonal block bit-for-bit (~1e-16/1e-17,
both director classes), confirming the extraction changed no numbers; (2) a
finite-difference check of addDirectorHessianRefNormalSens against
addDirectorHessianProduct's own value perturbed in t (~1e-11/1e-12, real
mode). Both formulas are closed-form polynomials in TacsScalar with no
fabs()/comparison/TacsRealPart()-gated branches, so they are holomorphic
and complex-step-safe by construction (SPEC's binding TACS_USE_COMPLEX
requirement) -- a full complex-mode libtacs.so rebuild was not run for this
specific test (real-mode regression-exact + FD checks were judged
sufficient given the manifestly polynomial, branch-free formulas; the a2d
op-level tests already exhaustively cover both real/complex modes for the
more novel second-order-A2D machinery).
…new ModRot test coverage; document a deeper-than-expected G1 obstruction

TACSDirector.h: TACSLinearizedRotation gets trivially-zeroing
addDirectorHessianProduct/addDirectorHessianRefNormalSens stubs (SPEC's own
"implementer's choice" allowance) so call sites needing to compile
generically across all three director classes don't need a C++17 if
constexpr (this codebase targets C++11).

tests/element_tests/shell_tests/test_beam_element.py: new
test_element_mat_sv_sens_nonlinear_director, bringing Beam2ModRot/
Beam3ModRot (TACSQuadraticRotation) into matrix-sens coverage via a
separate, narrow element list (not the shared self.elements, since
test_element_jacobian and others in this file are not yet analytic for
this director class and would break) -- SPEC-phase-7.md sec 6.3's own
flagged prerequisite for G1/G2's oracle to exist at all.

TACSBeamElement.h: attempted G1's TACS_STIFFNESS_MATRIX branch per
SPEC-phase-7.md sec 4.1's recipe (dd1phi/dd2phi/dd1psi/dd2psi built via
addDirectorHessianProduct, mirroring G2's own construction including its
tying-strain contribution) -- FAILED. The new test above caught two
independent, genuine problems, not implementation slips: (1) the recipe's
own "sweep rotational DOFs only" premise is incomplete -- a direct
re-derivation shows U contains a cubic (u, q, q)-permutation cross term
(u^T*Auq*B(q), B(q) quadratic in q via the director map) whose third
derivative is nonzero at TRANSLATIONAL output DOFs too, confirmed
empirically via a standalone getMatType-based diagnostic, independent of
TestElementMatSVSens's own FD; (2) the implemented ROTATIONAL-DOF term1+term2
computation has its own separate, unresolved mismatch. Documented in full at
the code site (TACSBeamElement.h) per this feature's attempt-first,
documented-failure protocol; TACS_STIFFNESS_MATRIX/TACS_MASS_MATRIX both
remain forwarded to the base FD/CS fallback for TACSQuadraticRotation/
TACSQuaternionRotation, unchanged in outward behavior from before this
task -- the new ModRot test coverage is currently trivially-passing FD-vs-FD
(honestly labeled as such) and becomes the regression net for a future
resolution.
…s concrete finding

Not attempted with new code this session -- SPEC-phase-7.md sec 3.2/4.2's
recipe ("contract the already-computed spsi/sphi stress-like weights,
used as the dd argument") is analogous in shape to G1's own recipe, which
Task 7.3 found was incomplete in a way only implementation-time
verification caught (a missing (u,q,q)-permutation cross term). Re-deriving
G2's actual need shows the missing "(dU/dC):(d^2C/dq^2)" term requires the
model's first-derivative weight evaluated at the REAL state (built from the
real vars-projected strain), not spsi/sphi (built from psi/phi-projected
kinematics) -- a genuinely new third kinematics chain
addMatXptSensInnerProductStiffness does not currently build at all, per its
own header comment. Given G1's directly-analogous "trust the dd-weight
framing literally" attempt produced a ~130%-off wrong result even after
adding the tying-strain correction, attempting G2's parallel construction
without first re-deriving and independently verifying the correct weight
(and checking for an equally hidden cross term) would very likely reproduce
the same failure mode. Documented at the code site with the full reasoning
trail; addMatXptSensInnerProduct's TACS_STIFFNESS_MATRIX branch is
unchanged, still forwarding TACSQuadraticRotation/TACSQuaternionRotation to
the base FD/CS fallback.
…sk 5.3's own finding from a new angle

G3 (addMatDVSensInnerProduct's TACS_GEOMETRIC_STIFFNESS_MATRIX branch): not
attempted with new code this session. Re-derived the "direct" piece from
scratch by tracing getMatType's GEOMETRIC branch precisely -- with
TACSLinearizedRotation's base strain identically zero, the Ju0x/Jd1x/Jd2x
Jacobian rows collapse to fixed constants and the resulting contraction
reduces exactly to epsi[0:4] (evalStrain's own axial/torsion/bending
components), making the "direct" piece's DV-sens a single
addGeometricTangentStressDVSens call. This CONFIRMS, rather than merely
repeats, Task 5.3's own finding: this is the same shape Task 5.3 already
tried and found insufficient, and the closure piece
(TacsBeamAddCrossDirectorJacobian/director::addDirectorJacobian) is a
genuinely separate "J(q)-congruence" pullback of a (u0xi,d01,d02)-level
Hessian block, distinct from the director-rate-level projection already
implicit in the direct piece. Given this session's own Task 7.3/7.4
experience (a recipe that looked complete on paper had a second,
independent missing term only caught by implementation-time verification),
attempting the closure piece without a full, independently re-verified
re-derivation was not started -- deferred with this added precision rather
than guessed at.

G4 (addMatXptSensInnerProduct's TACS_GEOMETRIC_STIFFNESS_MATRIX branch):
correspondingly blocked, per SPEC-phase-7.md sec 5's own hard sequencing
(G4 depends on G3 landing first).

No functional change -- both matType branches remain forwarded to the base
FD/CS fallback for all director classes, unchanged from Task 5.3/5.5's own
landed state.
@timryanb timryanb closed this Jul 29, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant