FESADev/docs/formulations/mitc4-shell-formulation.md

FESA Linear Static MITC4 Shell Formulation

Metadata

• feature_id: linear-static-mitc4-shell
• source_requirement: docs/requirements/linear-static-mitc4-shell.md
• source_research: docs/research/linear-static-mitc4-shell-research.md
• status: ready-for-numerical-review
• owner_agent: formulation-agent
• date: 2026-08-11
• product_execution_scope: small-strain, small-rotation linear static only
• future_formulation_scope: geometrically nonlinear Total Lagrangian residual/tangent; not executable

1. Purpose, authority, and decision boundary

This document is the mathematical and algorithmic contract for the FESA four-node MITC4 shell element. It strengthens the earlier draft with the approved requirements and research evidence. It does not define C++ APIs, parser classes, HDF5 paths, or a reference-solver workflow.

The normative product path is:

• one linear-static step;
• small strain and small rotation;
• four-node bilinear degenerated-continuum shell geometry;
• 20 physical element coordinates, embedded in 24 global shell DOFs;
• MITC4 edge-midpoint assumed transverse shear;
• homogeneous isotropic linear elasticity, one centered layer, constant thickness;
• one documented 2 x 2 x 2 integration path;
• nodal CLOAD and prescribed global shell DOFs only.

Section 15 records a future Total Lagrangian residual and consistent tangent. That section is formulation evidence only. It does not expand the current analysis procedure, mutable state, input subset, or implementation plan.

Source labels S4 and S4R both select this one FESA formulation by an approved semantic mapping. They do not select Abaqus integration or stabilization behavior, and this document makes no Abaqus formulation-equivalence claim.

Exact drilling scale/coefficient, drilling-energy warning ratio, smooth-director angle, drilling-load projection tolerance, and Jacobian/warpage thresholds remain Numerical Review decisions. Symbols are retained for them; no implementation default is implied.

2. Scope and assumptions

2.1 Included product formulation

• Reference coordinates: (xi, eta, zeta) in [-1,1]^3.
• Source node order: 1=(-1,-1), 2=(+1,-1), 3=(+1,+1), 4=(-1,+1).
• Global nodal DOF order: [UX, UY, UZ, URX, URY, URZ].
• Physical nodal fields: three midsurface translations and two director-tangent rotations.
• Numerical nodal field: one director-parallel drilling coordinate.
• Material: finite E > 0, -1 < nu < 0.5.
• Section: one centered homogeneous isotropic layer with finite constant t > 0.
• Plane-stress condition sigma_33 = 0; no thickness stretch.
• User-consistent units. Angles are dimensionless.
• Initial directors are derived from source connectivity and smooth-patch averaging; explicit normals are outside this feature.
• External loads are global nodal forces and non-drilling nodal moments.

2.2 Excluded product formulation

• Geometrically nonlinear execution, follower loads, dynamics, buckling, contact, material nonlinearity, and thermal strain.
• Composite or multilayer sections, offsets, variable/distributed thickness, orthotropy, and user material orientation.
• Pressure, gravity, body force, edge traction, DLOAD, or any distributed-load input.
• MITC4+, MITC4/D, EAS/ANS variants other than the shear tying defined here, reduced integration, hourglass control, and selective integration.
• A physical drilling strain, drilling stress/resultant, or drilling-direction applied nodal moment.
• Shared-node sharp folds, hinges, shell intersections, and shell-beam drilling transfer. A physical fold requires duplicate source nodes.

2.3 Kinematic and constitutive assumptions

The initial thickness director is a unit vector and thickness is a separate scalar. Straight material lines through the thickness remain straight in the physical five-DOF kinematics. The transverse normal strain and stress are excluded, while transverse shear deformation is retained. The shear correction factor is

kappa_s = \frac{5}{6}.

The 5/6 value is restricted to this homogeneous rectangular-thickness contract. It is not a general composite-shell rule.

3. Coordinates, signs, and location identity

3.1 Natural coordinates and positive face

The positive thickness direction is +zeta. With the approved node order, the positive element normal candidate follows

\mathbf n_e \parallel \mathbf A_\xi \times \mathbf A_\eta.

The bottom, middle, and top section positions are respectively

z=-t/2,\qquad z=0,\qquad z=+t/2,
\qquad z=\frac{t}{2}\zeta.

Positive section moment is defined by

\mathbf M=\int_{-t/2}^{t/2} z\,\boldsymbol\sigma_m(z)\,dz.

This convention fixes the curvature sign in Section 14. Reversing source connectivity reverses the positive face and is not silently repaired.

3.2 Integration-point local frame

At each midsurface location, let the interpolated positive-thickness director be

\overline{\mathbf d}=\frac{\sum_I N_I\mathbf d_I}
{\left\|\sum_I N_I\mathbf d_I\right\|}.

Then define

\mathbf e_3=\overline{\mathbf d},
\qquad
\mathbf e_1=
\frac{\mathbf A_\xi-(\mathbf A_\xi\cdot\mathbf e_3)\mathbf e_3}
     {\|\mathbf A_\xi-(\mathbf A_\xi\cdot\mathbf e_3)\mathbf e_3\|},

\mathbf e_2=\mathbf e_3\times\mathbf e_1.

Thus local direction 1 follows the projected +xi covariant tangent, direction 2 completes a right-handed frame, and direction 3 follows the interpolated positive thickness director. The source-order surface normal must have positive dot product with e_3; an opposed or degenerate pairing is invalid. The frame is evaluated and stored separately at every required shell location; one element-average frame is not substituted for warped geometry. Frame degeneracy is a geometry error.

All generalized strains, resultants, and section stresses in this document use this local frame. Nodal U/UR and RF/RM remain global.

3.3 Stable location identity

No recovery values at different locations are averaged to hide mismatches. A shell row is identified at least by source element identity, source element type, internal formulation identity, in-plane natural coordinates, section position when applicable, local frame/director, and component order. Exact HDF5 field names belong to the I/O contract.

4. Initial directors and nodal tangent frames

4.1 Element normal candidate

For element e, evaluate the midsurface covariant tangents at the element center:

\mathbf A_\xi^e(0,0)=\sum_{I=1}^4 N_{I,\xi}(0,0)\mathbf X_I,
\qquad
\mathbf A_\eta^e(0,0)=\sum_{I=1}^4 N_{I,\eta}(0,0)\mathbf X_I.

The source-order-positive candidate is

\mathbf n_e=
\frac{\mathbf A_\xi^e\times\mathbf A_\eta^e}
     {\|\mathbf A_\xi^e\times\mathbf A_\eta^e\|}.

Its deterministic surface-area weight is

A_e=\sum_{g=1}^{4}w_g
\left\|\mathbf A_\xi(\xi_g,\eta_g)
\times\mathbf A_\eta(\xi_g,\eta_g)\right\|,

using the same four midsurface Gauss locations as the element formulation.

4.2 Smooth shared-node director

For source node I and its incident accepted shell elements E(I), define

\mathbf s_I=\sum_{e\in E(I)} A_e\mathbf n_e,
\qquad
\mathbf d_I=\frac{\mathbf s_I}{\|\mathbf s_I\|}.

Incident elements are accumulated in stable source-element order. Before averaging, opposing orientations fail. After averaging, every incident deviation

\theta_{eI}=\cos^{-1}\!\left(
\operatorname{clamp}(\mathbf n_e\cdot\mathbf d_I,-1,1)\right)

must satisfy an approved smooth-patch bound theta_smooth. The research value 20 degrees is only the first Numerical Review candidate. It is not fixed here. Zero or near-zero cross products and averaged vectors fail using approved scale-aware tolerances. They are never replaced with a global axis.

4.3 Deterministic nodal tangent frame

Let (g_x,g_y,g_z) be the global Cartesian basis. Choose the least-aligned basis axis

\mathbf g_*=
\underset{\mathbf g\in\{\mathbf g_x,\mathbf g_y,\mathbf g_z\}}
{\arg\min}\ |\mathbf g\cdot\mathbf d_I|,

with ties resolved in the fixed order x, y, z. Then

\mathbf a_I=
\frac{\mathbf g_*-(\mathbf g_*\cdot\mathbf d_I)\mathbf d_I}
     {\|\mathbf g_*-(\mathbf g_*\cdot\mathbf d_I)\mathbf d_I\|},
\qquad
\mathbf b_I=\mathbf d_I\times\mathbf a_I.

The nodal frame (a_I,b_I,d_I) is orthonormal and right-handed: a_I x b_I=d_I.

5. Primary variables and six-to-five DOF embedding

5.1 Global and local nodal coordinates

At node I, the global shell vector is

\mathbf q_I^g=
\begin{bmatrix}
\mathbf u_I\\ \boldsymbol\theta_I^g
\end{bmatrix}
=
\begin{bmatrix}
U_X&U_Y&U_Z&UR_X&UR_Y&UR_Z
\end{bmatrix}^T.

Define

\mathbf R_I=
\begin{bmatrix}\mathbf a_I&\mathbf b_I&\mathbf d_I\end{bmatrix},
\qquad
\begin{bmatrix}\alpha_I\\\beta_I\\\gamma_I\end{bmatrix}
=\mathbf R_I^T\boldsymbol\theta_I^g.

alpha and beta are physical director-tangent rotations. gamma is the numerical drilling coordinate. Under the small-rotation convention

\delta\mathbf d_I
=\boldsymbol\theta_I^g\times\mathbf d_I
=\beta_I\mathbf a_I-\alpha_I\mathbf b_I,

and the gamma_I d_I part vanishes identically from director change.

5.2 Element transforms

Collect the four global nodal vectors as q_g in R^24. Define the block transform

\mathbf T=\operatorname{blockdiag}(\mathbf T_1,\ldots,\mathbf T_4),
\qquad
\mathbf T_I=
\begin{bmatrix}
\mathbf I_3&\mathbf 0\\
\mathbf 0&\mathbf R_I^T
\end{bmatrix}.

Let P in R^(20 x 24) select, for each node, [u_X,u_Y,u_Z,alpha,beta], and let D_gamma in R^(4 x 24) select the four gamma entries. Then

\boxed{\mathbf q_{20}=\mathbf T_p\mathbf q_g,
\qquad \mathbf T_p=\mathbf P\mathbf T},

\boxed{\boldsymbol\gamma=\mathbf T_d\mathbf q_g,
\qquad \mathbf T_d=\mathbf D_\gamma\mathbf T}.

T_p is 20 x 24 and T_d is 4 x 24. The physical strain operator acts only on q_20; gamma must not enter physical strain or recovery.

5.3 Virtual work preservation

For any physical local force f_20 and drill force m_d,

\delta\mathbf q_{20}^T\mathbf f_{20}
+\delta\boldsymbol\gamma^T\mathbf m_d
=\delta\mathbf q_g^T
\left(\mathbf T_p^T\mathbf f_{20}+\mathbf T_d^T\mathbf m_d\right).

This identity is the required transformation-energy check. A physical rigid mode is constructed from rigid midsurface translations and tangent-plane director changes; blindly assigning the complete spatial rigid-rotation vector, including its normal projection, to every UR entry would excite the numerical drill coordinate and is not the five-DOF rigid mode.

6. Strong form and boundary conditions

6.1 Three-dimensional parent statement

On the reference shell volume Omega_0, the linear elastic parent problem is

\nabla\cdot\boldsymbol\sigma=\mathbf0
\quad\text{in }\Omega_0,
\qquad
\boldsymbol\sigma=\mathbb C_{ps}:\boldsymbol\epsilon,

with

\boldsymbol\epsilon=\operatorname{sym}(\nabla\mathbf u),
\qquad \sigma_{33}=0.

The continuum strong form explains the element lineage. In the approved product path there is no body force and no surface/edge traction input.

6.2 Essential and natural data

Prescribed global shell DOFs form the Dirichlet data. The only external generalized forces are aggregated nodal CLOAD values. For a nodal moment M_I, define

M_{d,I}=\mathbf d_I\cdot\mathbf M_I.

If M_d,I is nonzero beyond an approved scale-aware projection tolerance, model validation fails as unsupported-drilling-load. The tolerance is a Numerical Review/I-O decision and is not replaced by an arbitrary absolute epsilon.

Distributed natural-boundary terms may appear in the parent continuum derivation, but their existence is not an executable load contract.

7. Weak form

7.1 Physical virtual work

Let V_20 be the admissible five-DOF element trial space and V_20^0 its homogeneous test space. The physical linear equilibrium is

\delta W_{int}^{phys}-\delta W_{ext}=0
\qquad\forall\delta\mathbf q_{20}\in V_{20}^0,

where

\delta W_{int}^{phys}
=\int_{\Omega_0}
\delta\overline{\boldsymbol\epsilon}^{,T}
\mathbf C_5
\overline{\boldsymbol\epsilon}\,dV,

and overbar denotes the MITC-projected shear field. In the current product path

\delta W_{ext}=\delta\mathbf q_g^T\mathbf f_e^{CLOAD}.

7.2 Numerical drilling potential

Drilling regularization is algebraically separate:

\Pi_d=\frac12\boldsymbol\gamma^T\mathbf K_d^l\boldsymbol\gamma,
\qquad
\delta\Pi_d=\delta\boldsymbol\gamma^T
\mathbf K_d^l\boldsymbol\gamma.

It does not approximate a shell strain energy and produces no physical generalized strain, resultant, or stress.

8. Bilinear discretization and degenerated geometry

8.1 Shape functions

With nodal signs (xi_I,eta_I),

N_I(\xi,\eta)=\frac14(1+\xi_I\xi)(1+\eta_I\eta),

or explicitly

\begin{aligned}
N_1&=\tfrac14(1-\xi)(1-\eta),&
N_2&=\tfrac14(1+\xi)(1-\eta),\\
N_3&=\tfrac14(1+\xi)(1+\eta),&
N_4&=\tfrac14(1-\xi)(1+\eta).
\end{aligned}

They must satisfy

\sum_I N_I=1,
\qquad N_I(\xi_J,\eta_J)=\delta_{IJ},
\qquad \sum_I N_{I,\xi}=\sum_I N_{I,\eta}=0.

8.2 Reference geometry

For element thickness t,

\boxed{
\mathbf X(\xi,\eta,\zeta)
=\sum_{I=1}^4N_I\mathbf X_I
+\frac{t\zeta}{2}\sum_{I=1}^4N_I\mathbf d_I
}.

The midsurface is X_0(xi,eta)=X(xi,eta,0). The initial covariant bases are

\mathbf G_\alpha=\frac{\partial\mathbf X}{\partial\xi^\alpha},
\qquad (\xi^1,\xi^2,\xi^3)=(\xi,\eta,\zeta).

8.3 Linear displacement interpolation

Using delta d_I=beta_I a_I-alpha_I b_I,

\boxed{
\mathbf u(\xi,\eta,\zeta)
=\sum_{I=1}^4N_I\mathbf u_I
+\frac{t\zeta}{2}\sum_{I=1}^4N_I
(\beta_I\mathbf a_I-\alpha_I\mathbf b_I)
}.

Therefore

\mathbf u=\mathbf H(\xi,\eta,\zeta)\mathbf q_{20},
\qquad \mathbf H\in\mathbb R^{3\times20}.

This interpolation reproduces constant translation. Its physical rigid-rotation test must combine u_I=omega x X_I with delta d_I=omega x d_I projected into (a_I,b_I). The corresponding global nodal rotation variable is the tangent projection

\boldsymbol\theta_I^g
=\boldsymbol\omega-(\boldsymbol\omega\cdot\mathbf d_I)\mathbf d_I,

so that gamma_I=0.

9. Mapping, Jacobian, and geometry validity

9.1 Covariant and contravariant bases

Define

\mathbf J=
\begin{bmatrix}\mathbf G_\xi&\mathbf G_\eta&\mathbf G_\zeta\end{bmatrix},
\qquad J=\det\mathbf J,

and contravariant bases by

\mathbf G^\alpha\cdot\mathbf G_\beta=\delta^\alpha_\beta.

The volume measure is

dV=J\,d\xi\,d\eta\,d\zeta.

9.2 Required validation locations

J and all basis quantities must be finite and valid at:

• all eight 2 x 2 x 2 stiffness points;
• the four midsurface MITC tying locations;
• any additional bottom/middle/top recovery evaluation location before recovery is committed.

Duplicate nodes, self-intersection, degenerate midsurface area, and reversed connectivity are separate fail-closed geometry errors.

9.3 Scale-aware measures

A dimensionless local determinant measure is

j_s=\frac{J}
{\|\mathbf G_\xi\|\,\|\mathbf G_\eta\|\,\|\mathbf G_\zeta\|}.

A valid point requires J>0 and j_s>tau_J; a valid element also satisfies an approved distortion/warpage inventory such as J_min/J_max. The exact tau_J and warpage bounds are Numerical Review decisions. No max(1, geometry_scale) clamp is permitted because it destroys consistent-unit scaling.

10. Linear kinematics and MITC4 shear projection

10.1 Direct covariant strain

For a physical element coordinate q_A, let

\mathbf h_A=\frac{\partial\mathbf u}{\partial q_A},
\qquad
\mathbf h_{A,\alpha}=\frac{\partial\mathbf h_A}{\partial\xi^\alpha}.

The direct covariant small-strain column is

B^{DI}_{\alpha\beta,A}
=\frac12\left(
\mathbf G_\alpha\cdot\mathbf h_{A,\beta}
+\mathbf G_\beta\cdot\mathbf h_{A,\alpha}
\right).

Equivalently,

\boldsymbol\epsilon^{DI}
=\operatorname{sym}\left(
\sum_{\alpha=1}^3
\mathbf u_{,\alpha}\otimes\mathbf G^\alpha
\right).

The transverse normal component is excluded by the shell constitutive contract.

10.2 Canonical tying locations

Use

T_{\xi-}=(0,-1,0),\qquad T_{\xi+}=(0,+1,0),

T_{\eta-}=(-1,0,0),\qquad T_{\eta+}=(+1,0,0).

The first pair supplies the covariant xi-zeta shear and the second pair supplies the covariant eta-zeta shear. Labels such as A/B/C/D are avoided because their assignment varies with published node conventions.

10.3 Assumed transverse shear

At any in-plane point,

\boxed{
\overline\epsilon_{\xi\zeta}(\xi,\eta)
=\frac{1-\eta}{2}\epsilon^{DI}_{\xi\zeta}(T_{\xi-})
+\frac{1+\eta}{2}\epsilon^{DI}_{\xi\zeta}(T_{\xi+})
},

\boxed{
\overline\epsilon_{\eta\zeta}(\xi,\eta)
=\frac{1-\xi}{2}\epsilon^{DI}_{\eta\zeta}(T_{\eta-})
+\frac{1+\xi}{2}\epsilon^{DI}_{\eta\zeta}(T_{\eta+})
}.

All in-plane covariant components remain direct. The assumed covariant tensor is reconstructed with the contravariant bases:

\overline{\boldsymbol\epsilon}
=\sum_{\alpha,\beta}
\overline\epsilon_{\alpha\beta}
\mathbf G^\alpha\otimes\mathbf G^\beta.

10.4 Local engineering strain vector

Transform the tensor to (e_1,e_2,e_3) and define

\boxed{
\overline{\mathbf e}=
\begin{bmatrix}
\epsilon_{11}&\epsilon_{22}&\gamma_{12}&\gamma_{13}&\gamma_{23}
\end{bmatrix}^T
=\overline{\mathbf B}\mathbf q_{20}},

where gamma_ij=2 epsilon_ij and B_bar in R^(5 x 20). The normative construction of B_bar is:

1. form each direct covariant strain column;
2. replace only its two transverse-shear covariant components by the same tying projection;
3. reconstruct the Cartesian tensor with the contravariant bases;
4. project it into the stored local frame and apply engineering-shear factors.

Strain, first variation, and future second variation must use the same linear MITC projection. Using tied shear in residual but direct shear in stiffness is forbidden.

11. Constitutive contract

11.1 Plane-stress material matrix

Let

G=\frac{E}{2(1+\nu)},
\qquad
\mathbf C_{ps}=\frac{E}{1-\nu^2}
\begin{bmatrix}
1&\nu&0\\
\nu&1&0\\
0&0&(1-\nu)/2
\end{bmatrix}.

For the engineering order of Section 10,

\boxed{
\mathbf C_5=
\begin{bmatrix}
\mathbf C_{ps}&\mathbf0\\
\mathbf0&\kappa_sG\mathbf I_2
\end{bmatrix}},

\mathbf s=
\begin{bmatrix}
\sigma_{11}&\sigma_{22}&\tau_{12}&\tau_{13}&\tau_{23}
\end{bmatrix}^T
=\mathbf C_5\overline{\mathbf e}.

C_5 is symmetric positive definite for the approved material range. sigma_33=0 is an assumption, not an emitted stress row.

There are no material state variables, history variables, or constitutive update in the approved linear elastic path.

11.2 Centered-section matrices as a cross-check

For an exactly linear in-plane strain through thickness,

\mathbf A=t\mathbf C_{ps},\qquad
\mathbf B=\mathbf0,\qquad
\mathbf D=\frac{t^3}{12}\mathbf C_{ps},\qquad
\mathbf A_s=\kappa_sGt\mathbf I_2.

These matrices are recovery and verification identities. The normative element stiffness remains the source-faithful volume integration in Section 12, so curved or warped mapping terms are not silently replaced by a flat analytical section kernel.

12. Linear element equations and drilling stabilization

12.1 Physical 20-DOF kernel

Using the quadrature of Section 13,

\boxed{
\mathbf K_{20}
=\int_{\Omega_0}\overline{\mathbf B}^T
\mathbf C_5\overline{\mathbf B}\,dV
\approx\sum_{g=1}^{8}
\overline{\mathbf B}_g^T\mathbf C_5\overline{\mathbf B}_g
J_gw_g},

\mathbf K_{20}\in\mathbb R^{20\times20},
\qquad
\mathbf f_{int}^{20}=\mathbf K_{20}\mathbf q_{20}.

The physical global contribution is

\boxed{
\mathbf K_{phys}^{24}=\mathbf T_p^T\mathbf K_{20}\mathbf T_p},
\qquad
\mathbf f_{phys}^{24}=\mathbf K_{phys}^{24}\mathbf q_g.

12.2 Drilling candidate contract

Let

\mathbf K_d^l=\operatorname{diag}(k_{d,1},k_{d,2},k_{d,3},k_{d,4}),
\qquad k_{d,I}>0,

with rotational-stiffness units force*length. Write

k_{d,I}=\alpha_d k_{ref,I},

where both the dimensionless coefficient alpha_d and the reference family remain unapproved until Numerical Review. Dimensionally compatible candidates carried from research are:

1. transverse-shear/area transition family

k_{ref,I}^{(A)}=
\frac{GtA_{eI}}{1+qA_{eI}/t^2},
\qquad q=2.5\times10^{-5},

where

A_{eI}=\int_{A_e}N_I\,dA
\approx\sum_{g=1}^{4}N_I(\xi_g,\eta_g)
\|\mathbf A_\xi\times\mathbf A_\eta\|_g w_g;

2. isotropic bending rigidity

k_{ref}^{(B)}=D_{iso}=\frac{Et^3}{12(1-\nu^2)};

3. a documented positive statistic formed only from the physical rotational block of K_20, whose entries all have force*length units.

The thesis rule 10^-3 min(all K_ii) is not admissible because it can mix translational force/length and rotational force*length diagonals. The order 10^-3 may only center a sensitivity study; it is not a nominal coefficient.

12.3 Stabilized 24-DOF matrix

The global drilling contribution and total element matrix are

\boxed{
\mathbf K_{drill}^{24}=\mathbf T_d^T\mathbf K_d^l\mathbf T_d},

\boxed{
\mathbf K_e^{24}=\mathbf K_{phys}^{24}+\mathbf K_{drill}^{24}},

\mathbf f_{int}^e=\mathbf K_e^{24}\mathbf q_g,
\qquad
\mathbf r_e=\mathbf f_{int}^e-\mathbf f_e^{CLOAD}.

No element distributed-load integral exists in the approved product routine. Mass and damping matrices are N/A for this linear-static feature.

For a valid free isolated element, the expected physical rank is 14. Embedding it in 24 coordinates creates the six physical rigid modes plus four drilling null modes. Four positive independent k_d,I values should remove only those drilling modes, giving expected stabilized rank 18 and nullity 6. These are verification targets, not substitutes for a computed singular-value/rank study.

12.4 Energy split

The element energies are

\boxed{
E_{phys}^e=\frac12\mathbf q_g^T
\mathbf K_{phys}^{24}\mathbf q_g
=\frac12\mathbf q_{20}^T\mathbf K_{20}\mathbf q_{20}},

\boxed{
E_{drill}^e=\frac12\mathbf q_g^T
\mathbf K_{drill}^{24}\mathbf q_g
=\frac12\boldsymbol\gamma^T\mathbf K_d^l\boldsymbol\gamma}.

Both have units force*length and are aggregated separately in stable source order. The ratio E_drill/E_phys is reported only when mathematically classifiable. If E_phys is zero or near zero, the two energies are reported explicitly; no arbitrary denominator clamp is used. The warning ratio remains open.

13. Numerical integration

13.1 Stiffness rule

Use tensor-product Gauss points

\xi,\eta,\zeta\in\left\{-\frac1{\sqrt3},+\frac1{\sqrt3}\right\},
\qquad w_\xi=w_\eta=w_\zeta=1.

Thus the physical kernel has four midsurface locations and two thickness locations, for eight volume evaluations. At each in-plane Gauss location:

1. evaluate and validate the midsurface local frame;
2. evaluate the four tying values once at zeta=0;
3. use the same tied covariant shear interpolation at both thickness points;
4. evaluate the remaining strain, mapping, material, and volume measure at the actual thickness point.

S4 and S4R use this identical rule. One-point reduced integration and Abaqus hourglass control are not selected by the S4R source label.

13.2 Prohibited quadrature substitutions

• replacing 2 x 2 midsurface integration with 1 x 1;
• using direct shear at any stiffness point;
• evaluating tying values from a different geometry or director inventory;
• averaging failed or missing tying values;
• changing the rule based on source S4 versus S4R;
• using higher-order quadrature in production without reopening Formulation and Numerical Review.

Higher-order or analytical integration may be used only as an independent verification cross-check.

14. Output recovery

14.1 Nodal solution, reaction, and equilibrium

Store at every source node, in global coordinates,

[U1,U2,U3,UR1,UR2,UR3]

and

[RF1,RF2,RF3,RM1,RM2,RM3].

The reaction/equilibrium vector is the assembled full residual

\mathbf r=\mathbf K\mathbf d-\mathbf F.

Constrained components are physical reactions; free components are residual evidence. Reaction is not reconstructed by summing section resultants.

For a fixed global reference point X_o, the force and moment balance evidence is

\mathbf e_F=\sum_I\left(\mathbf F_I^{CLOAD}+\mathbf R_I\right),

\mathbf e_M=\sum_I\left[
(\mathbf X_I-\mathbf X_o)\times
(\mathbf F_I^{CLOAD}+\mathbf R_I)
+\mathbf M_I^{CLOAD}+\mathbf{RM}_I\right].

Only constrained residual components contribute to R/RM; free components remain separate residual evidence. Normalization and acceptance thresholds belong to Numerical Review and the I/O verification contract.

14.2 Generalized strain by thickness moments

At each of the four midsurface Gauss locations, let

\mathbf e_m(z)=
\begin{bmatrix}\epsilon_{11}&\epsilon_{22}&\gamma_{12}\end{bmatrix}^T,
\qquad
\boldsymbol\gamma_s(z)=
\begin{bmatrix}\gamma_{13}&\gamma_{23}\end{bmatrix}^T,

evaluated in the same stored local frame. Define

\boldsymbol\epsilon_0=
\frac1t\int_{-t/2}^{t/2}\mathbf e_m(z)\,dz,

\boldsymbol\kappa=
\frac{12}{t^3}\int_{-t/2}^{t/2}z\,\mathbf e_m(z)\,dz,

\boldsymbol\gamma_0=
\frac1t\int_{-t/2}^{t/2}\boldsymbol\gamma_s(z)\,dz.

Use the same two-point thickness quadrature:

\int_{-t/2}^{t/2}f(z)\,dz
\approx\frac t2\sum_{h=1}^2w_h f(t\zeta_h/2).

The mandatory generalized strain order is

\boxed{
[E11,E22,G12,K11,K22,K12,G13,G23]
=[\boldsymbol\epsilon_0,\boldsymbol\kappa,
\boldsymbol\gamma_0]}.

Membrane and shear strains are dimensionless; curvature has units 1/length. For an exactly linear in-plane field, e_m(z)=epsilon_0+z kappa, fixing the sign convention.

14.3 Section resultants

With

\boldsymbol\sigma_m=
\begin{bmatrix}\sigma_{11}&\sigma_{22}&\tau_{12}\end{bmatrix}^T,
\qquad
\boldsymbol\tau_s=
\begin{bmatrix}\tau_{13}&\tau_{23}\end{bmatrix}^T,

define

\mathbf N=\int_{-t/2}^{t/2}\boldsymbol\sigma_m\,dz,
\qquad
\mathbf M=\int_{-t/2}^{t/2}z\boldsymbol\sigma_m\,dz,
\qquad
\mathbf Q=\int_{-t/2}^{t/2}\boldsymbol\tau_s\,dz.

The mandatory order is

\boxed{
[N11,N22,N12,M11,M22,M12,Q13,Q23]
=[\mathbf N,\mathbf M,\mathbf Q]}.

N and Q have units force/length; M has units force. For the centered homogeneous linear-through-thickness case,

\mathbf N=\mathbf A\boldsymbol\epsilon_0,
\qquad
\mathbf M=\mathbf D\boldsymbol\kappa,
\qquad
\mathbf Q=\mathbf A_s\boldsymbol\gamma_0,

which is a required recovery cross-check.

14.4 Bottom, middle, and top in-plane stress

At each midsurface Gauss location, perform additional recovery evaluations at zeta=-1,0,+1 and emit only

\boxed{[S11,S22,S12]}.

These are direct section-position evaluations, not extrapolations from a different location. S33=0 is documented but not emitted. S13/S23 point stresses are not emitted because the approved external contract carries transverse shear through Q13/Q23.

14.5 No nodal extrapolation

Element generalized quantities and section stresses remain at their defined shell locations. This formulation does not define nodal averaging or extrapolation.

15. Future-only geometrically nonlinear residual and tangent

Non-executable boundary: This section records the requested future formulation. It must not be routed through the current linear-static procedure or state.

15.1 Total Lagrangian kinematics

Let the physical nonlinear coordinate vector retain three translations and two director parameters per node. In the reference configuration,

\mathbf x(\xi,\eta,\zeta;\mathbf q)
=\sum_I N_I(\mathbf X_I+\mathbf u_I)
+\frac{t\zeta}{2}\sum_I N_I\mathbf d_I(\boldsymbol\phi_I),

where d_I(phi_I) is a unit-director update defined by an approved finite-rotation map. A candidate two-parameter exponential chart is

\mathbf d_I(\boldsymbol\phi_I)=
\exp\!\left([
\phi_{1I}\mathbf a_I+\phi_{2I}\mathbf b_I]_{\times}\right)
\mathbf d_I^0.

The deformation gradient and Green-Lagrange strain are

\mathbf F=\mathbf x_{,\alpha}\otimes\mathbf G^\alpha,
\qquad
\mathbf E=\frac12(\mathbf F^T\mathbf F-\mathbf I).

For physical coordinates q_A and q_B, define

\mathbf h_A=\frac{\partial\mathbf x}{\partial q_A},
\qquad
\mathbf h_{AB}=\frac{\partial^2\mathbf x}
{\partial q_A\partial q_B}.

The covariant Green-Lagrange components and their exact first and second derivatives are

E_{\alpha\beta}=\frac12\left(
\mathbf x_{,\alpha}\cdot\mathbf x_{,\beta}
-\mathbf G_\alpha\cdot\mathbf G_\beta\right),

\boxed{
E_{\alpha\beta,A}=\frac12\left(
\mathbf h_{A,\alpha}\cdot\mathbf x_{,\beta}
+\mathbf x_{,\alpha}\cdot\mathbf h_{A,\beta}\right)},

\boxed{
\begin{aligned}
E_{\alpha\beta,AB}=\frac12\big(&
\mathbf h_{AB,\alpha}\cdot\mathbf x_{,\beta}
+\mathbf x_{,\alpha}\cdot\mathbf h_{AB,\beta}\\
&+\mathbf h_{A,\alpha}\cdot\mathbf h_{B,\beta}
+\mathbf h_{B,\alpha}\cdot\mathbf h_{A,\beta}
\big).
\end{aligned}}

For translational coordinates, h_AB=0. Rotation-rotation blocks retain the exact second derivative of the director map; mixed-node blocks of h_AB are zero for the nodal interpolation above. The quadratic products of h_A and h_B remain and form the stress-dependent geometric contribution even when h_AB=0.

Apply the same covariant MITC projection to the two transverse shear components of E, its first variation, and its second variation:

\overline{\mathbf E}=\mathcal P_{MITC}(\mathbf E),

\delta\overline{\mathbf E}
=\mathcal P_{MITC}(\delta\mathbf E),
\qquad
\delta\Delta\overline{\mathbf E}
=\mathcal P_{MITC}(\delta\Delta\mathbf E).

The exact first and second derivatives of the chosen director map are mandatory. The mixed director Hessian is not assumed zero merely because the linearized small-rotation expression is linear.

Tensor-to-engineering transformation for the future material law uses the fixed initial local frame at each reference integration point. Consequently B_A and G_AB below are obtained by applying, in order, the MITC projection, the fixed reference-frame tensor transformation, and the engineering-shear factors to the explicit derivatives above. A future rotating output frame is a recovery decision and must not be inserted into the constitutive tangent without its own derivatives.

15.2 Residual

Use second Piola-Kirchhoff stress S work-conjugate to E. In engineering-vector notation, define

\mathbf B_A(\mathbf q)=
\frac{\partial\overline{\mathbf e}}{\partial q_A}.

For the geometrically nonlinear but materially linear candidate,

\mathbf s=\mathbf C_5\overline{\mathbf e},
\qquad
\mathbf C_T=\mathbf C_5,

interpreted as a second-Piola/Green-Lagrange St. Venant-Kirchhoff-type shell law in the fixed reference frame with the same plane-stress and shear-correction boundary as the linear kernel. A different finite-strain or material-nonlinear law requires a new constitutive formulation decision.

Then the physical internal residual is

\boxed{
r_A^{phys}(\mathbf q)
=\int_{\Omega_0}\mathbf B_A^T\mathbf s\,dV}.

The total residual is

\mathbf r(\mathbf q)=
\mathbf r^{phys}(\mathbf q)
+\mathbf r^{drill}(\mathbf q)-\mathbf f_{ext}(\mathbf q).

15.3 Consistent material and geometric tangent

Let

\mathbf C_T=\frac{\partial\mathbf s}
{\partial\overline{\mathbf e}},
\qquad
\mathbf G_{AB}=\frac{\partial^2\overline{\mathbf e}}
{\partial q_A\partial q_B}.

The consistent physical tangent is

\boxed{
K_{AB}^{phys}=K_{AB}^{mat}+K_{AB}^{geo}},

\boxed{
K_{AB}^{mat}=\int_{\Omega_0}
\mathbf B_A^T\mathbf C_T\mathbf B_B\,dV},

\boxed{
K_{AB}^{geo}=\int_{\Omega_0}
\mathbf s^T\mathbf G_{AB}\,dV}.

K_geo contains the stress-dependent second variation of Green-Lagrange strain, including consistent finite-director derivatives. The future total tangent is

\boxed{
\mathbf K_T=
\mathbf K^{mat}+\mathbf K^{geo}+\mathbf K^{drill}-\mathbf K^{load}}.

For a constant generalized dead-load vector in the selected coordinate chart, K_load=0. Follower pressure and other configuration-dependent loads are outside scope. Finite-rotation nodal-moment work must be separately defined before assuming its load tangent is zero.

15.4 Future drilling tangent boundary

If the linear initial-frame penalty is formally retained,

\Pi_d=\frac12\boldsymbol\gamma(\mathbf q)^T
\mathbf K_d^l(\mathbf q)\boldsymbol\gamma(\mathbf q),
\qquad
\mathbf r^{drill}=\frac{\partial\Pi_d}{\partial\mathbf q},
\qquad
\mathbf K^{drill}=\frac{\partial^2\Pi_d}{\partial\mathbf q^2}.

However, a constant initial-frame diagonal penalty is not automatically objective under large rotation. The nonlinear frame update, drilling coordinate, scale update, and their first/second derivatives are unresolved future decisions. They do not block the current linear kernel but do block geometrically nonlinear implementation.

15.5 Newton equation and linear limit

At iteration k,

\mathbf K_T(\mathbf q^{(k)})\Delta\mathbf q
=-\mathbf r(\mathbf q^{(k)}),
\qquad
\mathbf q^{(k+1)}=\mathbf q^{(k)}+\Delta\mathbf q.

At the undeformed, stress-free state with the small-rotation director derivative,

\mathbf K^{geo}=\mathbf0,

and the physical material tangent reduces to K_phys^24 in Section 12. With the approved linear drill block, the complete limit is K_e^24.

16. Algorithm pseudocode

16.1 Deterministic preprocessing

for each shell element in stable source order:
    validate four distinct source nodes and source-order geometry
    compute center positive normal candidate n_e
    compute 2x2 surface area A_e

for each shell source node in stable source order:
    gather incident candidates in stable element order
    reject opposed, degenerate, or too-sharp incident normals
    d_I = normalize(sum(A_e * n_e))
    select least-aligned global axis with deterministic tie break
    construct right-handed (a_I, b_I, d_I)

16.2 Linear element kernel

input: X_I, d_I, a_I, b_I, E, nu, t, q_g, nodal CLOAD share
build T, T_p, T_d
initialize K20[20,20] = 0

evaluate and validate four midsurface tying locations
for each 2x2 midsurface Gauss location in fixed order:
    construct and validate local frame (e1,e2,e3)
    compute the four direct covariant tying shear B columns at zeta=0
    for each thickness Gauss point in fixed order:
        build degenerated geometry, J, reciprocal bases, H derivatives
        validate finite positive scale-aware Jacobian
        form direct covariant B columns
        replace only xi-zeta and eta-zeta shear by MITC tying interpolation
        transform to local engineering B_bar[5,20]
        K20 += B_bar^T * C5 * B_bar * J * weight

check K20 finite and symmetric within approved normalized tolerance
Kphys24 = T_p^T * K20 * T_p
construct candidate-aware positive Kd_local without closing open coefficient choices
Kdrill24 = T_d^T * Kd_local * T_d
Ke24 = Kphys24 + Kdrill24
fint24 = Ke24 * q_g
residual24 = fint24 - f_CLOAD
return matrices, residual, transforms, frames, and separate energy operators

16.3 Global linear-static lifecycle

assemble all Ke24 contributions with stable element-local COO ordering
partition full K into Kff, Kfc, Kcf, Kcc in stable free/constrained order
factorize Kff before load assembly
assemble and deterministically aggregate nodal CLOAD
solve Kff * df = Ff - Kfc * dc
reconstruct full displacement d
compute full residual r = K*d - F
recover shell rows and physical/drilling energies in stable source order
validate complete finite candidate state/output, then commit

16.4 Recovery

for each element and each 2x2 midsurface location in fixed order:
    evaluate two thickness Gauss strains/stresses with the stiffness formulation
    integrate thickness zeroth/first moments for generalized strain/resultants
    independently evaluate S11,S22,S12 at zeta = -1, 0, +1
    attach exact natural coordinates, section position, frame, and source identity

recover nodal global U/UR and full-residual RF/RM
compute E_physical and E_drill separately; never clamp a near-zero denominator

16.5 Future nonlinear tangent check

given a future approved finite-director map and nonlinear drill potential:
    evaluate residual r(q)
    evaluate analytic consistent K_T(q)
    for several normalized perturbation directions p and decreasing h:
        compare K_T*p with [r(q+h*p)-r(q-h*p)]/(2*h)
    require the error to decrease in the expected truncation range

17. Verification contract

17.1 Algebraic and geometry checks

• Shape-function partition of unity, Kronecker delta, and derivative sums.
• Nodal and integration frames orthonormal and right-handed.
• Finite positive Jacobian at every Gauss, tying, and committed recovery point.
• K20, Kphys24, Kdrill24, and Ke24 symmetry.
• Transformation work/energy invariance.
• Physical rigid modes satisfy normalized Kphys24 r=0.
• Stabilized free-element nullity is exactly six; accepted non-rigid physical modes have positive physical energy.
• Pure drill vectors have zero physical energy and positive drilling energy.
• Consistent force/length unit rescaling leaves dimensionless decisions unchanged.

Requirements currently prescribe normalized targets 1e-12 for frame, symmetry, and transformation-energy checks and 1e-10 for rigid action, residual, and global equilibrium unless Numerical Review documents replacements.

17.2 Patch and sign checks

Independently verify:

• constant membrane strain/stress and N sign;
• pure bending about both local axes, K11/K22, M, and bottom/top stress sign;
• constant transverse shear and Q13/Q23 order;
• pure twist and K12/M12 convention;
• zero strain/resultant/stress contribution from a pure drilling vector.

17.3 Locking, distortion, and curved shells

• Thin and thick plate/shell mesh and thickness sequences are required; one displacement on one mesh is insufficient.
• Distorted and warped valid quadrilaterals must be swept through approved geometry measures.
• Original MITC4 controls transverse-shear locking but can retain membrane locking in distorted curved meshes. This is a known limitation, not permission to add MITC4+.
• Preferred nodal-load-compatible curved benchmarks are the pinched cylinder and NAFEMS LE3 hemispherical shell. Scordelis-Lo is admissible only after an equivalent nodal-load adaptation is documented.

17.4 Drilling sensitivity

For every candidate reference scale, sweep dimensionless coefficients around the research order of magnitude and record:

• free-element rank and constrained system conditioning;
• global U and physical N/M/Q sensitivity;
• E_phys and E_drill without denominator clamping;
• invariance under consistent unit conversion, thickness ratios, and mesh sizes.

The study must select the reference family, nominal coefficient/range, and warning criterion before Implementation Planning.

17.5 Reference-comparison boundary

Abaqus comparisons block only on matched global U1/U2/U3 rows under the approved mixed tolerance decided downstream. UR1/UR2/UR3 is fully reported and may emit a deterministic nonblocking large-error warning, but it does not change pass/fail. FESA S4 and S4R inputs must produce the same internal numerical rows for identical supported models while preserving distinct source metadata. Abaqus S4 and S4R are not expected to be numerically identical on finite meshes.

17.6 Future nonlinear tangent verification

Before any geometrically nonlinear implementation, verify residual directional derivatives, tangent symmetry for conservative loading, objectivity under large rigid motion, zero stress-free geometric stiffness, and convergence of Newton iterations. These checks are future-only and do not authorize a nonlinear procedure.

18. Numerical risks

Risk	Consequence	Required control
Transverse-shear locking	overly stiff thin-shell response	exact edge-midpoint MITC projection and thickness/mesh convergence
Membrane locking on distorted curved meshes	slow or nonuniform convergence	distortion/curvature sweeps; document original MITC4 limitation
Volumetric locking	N/A for the approved plane-stress shell contract	do not infer a three-dimensional incompressible formulation
Wrong tying pair or engineering-shear factor	swapped/incorrect shear and loss of patch consistency	component-level tying and patch tests
Reversed or degenerate Jacobian	invalid basis, sign, or energy	fail at every required point using scale-aware measures
Discontinuous shared director	artificial coupling or undefined frame	reject; require duplicate nodes at folds
Fixed-axis tangent singularity	nondeterministic rotation transform	least-aligned global-axis construction
Drilling coefficient too small	rank/conditioning failure	rank and conditioning sensitivity sweep
Drilling coefficient too large	contaminated displacement/resultant	physical-output and separate-energy sensitivity sweep
Mixed-unit diagonal scale	unit-dependent stiffness	prohibit min(all Kii); use rotationally compatible scale
Misconstructed rigid test	false drill energy in a physical mode	use rigid translation plus tangent director change with gamma=0
Recovery/stiffness mismatch	inconsistent energy and section output	same frames, tying, material, and thickness rule
Location averaging	hidden sign/identity error	preserve exact location rows; no nodal extrapolation
Treating S4R as reduced FESA integration	divergent kernel and unsupported hourglass behavior	one documented rule for S4/S4R
Nonlinear director Hessian omission	inconsistent future tangent	exact first/second derivative of approved rotation map
Nonobjective nonlinear drill penalty	artificial large-rotation energy	future review of rotating frame and drill potential

19. Evidence basis and applicability

19.1 Repository research sources

The detailed source tiers, extracted facts, page references, benchmark provenance, and evidence limits are owned by docs/research/linear-static-mitc4-shell-research.md. The primary local source set under docs/reference-papers/MITC4/ includes:

• AContinuumMechanicsBasedFourNodeShell_001.md and _002.md;
• FourNodeQuadrilateralShellElementMITC4_001.md;
• MITC공부_001.md and _002.md;
• 유한요소해석법을이용한쉘구조물의동적좌굴해석_001.md through _008.md;
• 쉘구조물의유한요소해석에대하여_001.md and _002.md.

The evidence-backed core is the degenerated-continuum geometry, 20 physical DOFs, edge-midpoint covariant shear tying, plane-stress constitutive response, and 2 x 2 midsurface plus two-point thickness integration. Six-global-DOF drilling regularization is a FESA interface/stability layer and must not be described as the physical MITC4 strain field.

19.2 FEM wiki cross-checks

• [[MITC4 Shell Element]]: four-node degenerated-continuum kinematics, two director-tangent rotations, and edge-midpoint assumed shear.
• [[MITC Shell Kinematics]]: midsurface/director/thickness interpolation.
• [[Assumed Transverse Shear Strain Interpolation]]: locking mechanism and tying rationale.
• [[Continuum Mechanics Based Four-Node Shell Element]]: three-dimensional virtual work lineage and benchmark thread.
• [[Shell Locking Phenomenon]]: separate transverse-shear and curved-shell membrane locking risks.
• [[Total Lagrangian Shell Formulation]] and [[Green-Lagrange Strain Linearization]]: future residual and material/geometric tangent separation.

Wiki pages are navigation and cross-check evidence; the approved requirements and research brief remain the project source of truth.

20. Requirement traceability

Requirement group	Formulation coverage	Remaining owner
001-004, 030, 037	linear-static boundary; S4/S4R one FESA path, source identity distinct	I/O, planning
005, 031-038	24 global DOFs; 20 physical plus four drill coordinates; transform, rank, energy	Numerical Review for drill choices
006-010	isotropic plane stress, one centered constant-thickness layer	I/O validation
011-016	connectivity-positive normals, deterministic area weighting, frames, Jacobian checks	Numerical Review for thresholds
017-020	global BC/CLOAD virtual work; parallel moment rejected; distributed loads excluded	I/O projection diagnostic
024-029	deterministic element buffers, partitioned linear lifecycle, full-residual reaction	planning
039-048	nodal/global and shell/local recovery inventory, units, identities, energy split	I/O schema
049-057	invariants, patches, locking, geometry, curved shells, drilling, equilibrium	Numerical Review/reference/physics
058-064	U blocking, UR warning-only comparison boundary	Numerical Review/reference model
065-072	reference immutability and displacement-only evidence boundary acknowledged	reference model

21. Open issues and downstream handoff

21.1 Blocking Numerical Review decisions

1. Select the dimensionally compatible drilling reference-scale family.
2. Select alpha_d, its sensitivity range, conditioning/rank acceptance, and the physical-output contamination bound.
3. Define classification and warning behavior for E_drill/E_phys, including the zero/near-zero physical-energy case.
4. Approve theta_smooth after curved-mesh resolution sweeps; 20 degrees is only the initial candidate.
5. Approve dimensionless Jacobian, distortion, and warpage measures and thresholds.
6. Approve the scale-aware tolerance used to reject a director-parallel nodal moment.
7. Re-derive/check the exact B_bar component signs, tying-point inventory, engineering-shear factors, rank, and bottom/top stress sign independently.

Until these items are approved, the formulation is ready for Numerical Review but not for Implementation Planning.

21.2 I/O Definition handoff

• Preserve source S4/S4R separately from internal FESA-MITC4.
• Define exact keyword subset, section/material resolution, and fail-closed diagnostics for director, folds, Jacobians, unsupported loads, and recovery.
• Define exact HDF5 row schemas for global U/UR, RF/RM; four midsurface generalized-strain/resultant locations; bottom/middle/top stress positions; full residual/equilibrium; and separate energies.
• Preserve local frame and natural-coordinate identity without location averaging.

21.3 Reference Model handoff

• Use at least one S4 and one S4R source artifact, but compare formulation-independent global displacement evidence rather than claiming element equivalence.
• Prioritize pinched cylinder and NAFEMS LE3 models compatible with nodal loads and approved BC semantics.
• Propose the mixed U tolerance and nonblocking UR warning threshold.
• Do not create, repair, rename, or run reference artifacts during this formulation gate.

21.4 Implementation Planning handoff

• Do not begin until Numerical Review closes Section 21.1.
• Translate the deterministic preprocessing, 24-to-20 transform, tying projection, physical/drill split, quadrature, recovery, and invariant portfolio into RED -> GREEN -> VERIFY tests before production changes.
• Keep future nonlinear state and tangent out of the current linear-static plan.

21.5 Future nonlinear formulation handoff

Before geometric-nonlinear implementation, separately approve the finite-director chart and update law, objective drilling potential and scale update, exact first/second director derivatives, nonlinear constitutive boundary, nodal-moment work/load tangent, state ownership, convergence controls, and nonlinear output contract. Section 15 supplies the residual/tangent structure but deliberately does not close those product decisions.