Skip to content

fix: use consistent Voigt order for the shell reference curvature - #97

Merged
oberbichler merged 1 commit into
masterfrom
fix/shell-curvature-order
Jul 26, 2026
Merged

fix: use consistent Voigt order for the shell reference curvature#97
oberbichler merged 1 commit into
masterfrom
fix/shell-curvature-order

Conversation

@oberbichler

Copy link
Copy Markdown
Owner

Problem

IgaShell3PAD builds the membrane metric (ref_a/act_a) and the strain transformation matrix T in Voigt order [11, 22, 12], and the actual curvature uses the same order:

act_b(act_a1_1.dot(act_a3), act_a2_2.dot(act_a3), act_a1_2.dot(act_a3));   // [11, 22, 12]

But the reference curvature built in add() used [11, 12, 22]:

ref_b(ref_a1_1.dot(ref_a3), ref_a1_2.dot(ref_a3), ref_a2_2.dot(ref_a3));   // [11, 12, 22]  <-- swapped

kap = T * (act_b - ref_b) applies the same [11, 22, 12] transform to both operands, so the mismatched ordering subtracts the b12 and b22 components across different slots — a wrong bending strain whenever the reference surface is curved (ref_b != 0).

Flat references (ref_b == 0) and the membrane element (IgaMembrane3PAD, whose ref_a/act_a are consistently [11, 22, 12]) are unaffected. This is why the existing baked reference test — a flat patch — does not change.

Fix

Reorder ref_b to [11, 22, 12] to match act_b and the transformation:

ref_b(ref_a1_1.dot(ref_a3), ref_a2_2.dot(ref_a3), ref_a1_2.dot(ref_a3));

Why the usual tests can't catch this

  • The baked test_iga_shell_3p_ad asserts values produced by the element itself, and its reference is flat — the bug is invisible there.
  • A finite-difference check of g vs f cannot catch it either: the AD differentiates whatever energy expression is written, so g stays consistent with f even when that expression is wrong.

Test

tests/test_iga_shell_3p_symmetry.py uses a value-level physical invariant on a curved reference: for an isotropic material the shell energy must be invariant under relabelling the parametric directions u1 <-> u2 (which only permutes the shape-function rows).

  • Before the fix: f_orig = 0.4182 vs f_swapped = 0.3031 → fails.
  • After the fix: equal → passes.

Cross-checked independently against a numpy reference implementation using the consistent [11, 22, 12] ordering (element f moves from 0.4182 (matches the buggy ordering) to 0.2237 (matches the reference)). Full suite: 44 passed, existing shell baseline unchanged.

Impact

Kirchhoff-Love shells with a curved reference geometry previously produced an incorrect energy, gradient and Hessian. Flat references and membranes were not affected.

IgaShell3PAD builds the membrane metric (ref_a/act_a) and the strain
transformation matrix in Voigt order [11, 22, 12], and the actual
curvature act_b uses the same order:

    act_b = (a1_1 . a3, a2_2 . a3, a1_2 . a3)

The reference curvature computed in add(), however, used [11, 12, 22]:

    ref_b = (a1_1 . a3, a1_2 . a3, a2_2 . a3)

Since kap = T * (act_b - ref_b) applies the same [11, 22, 12] transform
to both, the mismatched ordering mixes the b22 and b12 curvature
components, giving a wrong bending strain whenever the reference surface
is curved (ref_b != 0). Flat references (ref_b == 0) and the membrane
element are unaffected, which is why the baked reference test (a flat
patch) does not change.

Reorder ref_b to [11, 22, 12] to match act_b and the transformation.

Add a value-level regression test: on a curved reference the isotropic
shell energy must be invariant under relabelling the parametric
directions (u1 <-> u2, i.e. permuting the shape-function rows). This
fails before the fix and passes after; it is verified independently
against a numpy reference implementation using the consistent ordering.
@oberbichler
oberbichler merged commit 75aefd7 into master Jul 26, 2026
10 checks passed
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