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`
- source_numerical_review: `docs/numerical-reviews/linear-static-mitc4-shell-review.md`
- status: `ready-for-numerical-review`
- owner_agent: `formulation-agent`
- date: `2026-08-12`
- revision_basis: approved independent-reference policy, fixed drilling rule, and
  removal of calibration gates `NR-O01` through `NR-O04`
- revision_state: `ready-for-numerical-rereview-not-implementation-planning`
- 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.

The drilling rule is fixed below by the approved requirements. Coefficient sweeps,
drilling-energy warnings and datasets, smooth-director calibration (`NR-O03`), and
distortion/warp threshold calibration (`NR-O04`) are outside the implementation gate.
The drilling-load projection and normalized algebraic checks remain fixed below.

## 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`, sort its incident accepted shell elements `E(I)` by stable
source-element identity. Before any averaging, require every unordered incident
pair to have the same orientation hemisphere:

$$
c_{pair,I}=\min_{e<f,\ e,f\in E(I)}\mathbf n_e\cdot\mathbf n_f>0.
$$

For a node with only one incident element this pair test is vacuous. A nonpositive
pair dot product is an `opposed-incident-normal` geometry error; it is not repaired
by flipping a candidate or by letting cancellation occur in the average. This
pairwise rule is deterministic and stronger than testing only the final average.

After the pairwise test, 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 the same stable order. No additional
`theta_smooth` rejection is applied. Nonfinite or zero cross products and averaged
vectors fail; they are never replaced with a global axis. The pairwise orientation
rule above remains the exact supported-patch predicate.

### 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

$$
\rho_{M,I}=\frac{|\mathbf d_I\cdot\mathbf M_I|}{\|\mathbf M_I\|}.
$$

An exactly zero aggregated moment is accepted as a separate zero case and does not
form this ratio. Every nonzero moment is accepted only when

$$
\boxed{\rho_{M,I}\le 10^{-12}}.
$$

Otherwise model validation fails as `unsupported-drilling-load`. The denominator is
never clamped. This is a direction/projection check and is invariant under a
consistent moment-unit conversion.

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_24` be the admissible global six-DOF element trial space after essential
boundary data are applied, and let `V_24^0` be its homogeneous test space. With

$$
\mathbf q_{20}=\mathbf T_p\mathbf q_g,
\qquad
\delta\mathbf q_{20}=\mathbf T_p\delta\mathbf q_g,
$$

the MITC-projected physical internal virtual work is

$$
\delta W_{int}^{phys}
=\int_{\Omega_0}
(\delta\overline{\boldsymbol\epsilon})^T
\mathbf C_5
\overline{\boldsymbol\epsilon}\,dV
=\delta\mathbf q_{20}^T\mathbf f_{int}^{20}
=\delta\mathbf q_g^T\mathbf T_p^T\mathbf f_{int}^{20},
$$

where

$$
\mathbf f_{int}^{20}
=\int_{\Omega_0}\overline{\mathbf B}^T
\mathbf C_5\overline{\boldsymbol\epsilon}\,dV.
$$

The overbar denotes the common MITC projection used by strain, residual, and
stiffness. The current product external virtual work is

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

### 7.2 Numerical drilling potential and complete weak form

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.

The complete linear 24-DOF equilibrium statement is therefore

$$
\boxed{
\delta\mathbf q_g^T\left[
\mathbf T_p^T\mathbf f_{int}^{20}
+\mathbf T_d^T\mathbf K_d^l\boldsymbol\gamma
-\mathbf f_e^{CLOAD}
\right]=0
\qquad\forall\delta\mathbf q_g\in V_{24}^0.}
$$

This is the normative weak form. A physical-only statement in `V_20` is merely the
restriction of this equation and must not be combined with a 24-DOF external-work
term as if both variations occupied the same test space.

## 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.

The center `(0,0,0)` is also included in the orientation and warpage inventory.
Duplicate nodes, self-intersection, degenerate midsurface area, and reversed
connectivity are separate fail-closed geometry errors. No failed location is
discarded or replaced by a value from another point.

### 9.3 Basic geometry predicates

Let the consecutive midsurface edge inventory be

$$
\mathcal E=\{(1,2),(2,3),(3,4),(4,1)\},
\qquad
\boxed{L_e=\max_{(I,J)\in\mathcal E}\|\mathbf X_J-\mathbf X_I\|}.
$$

`L_e` must be finite and strictly positive. At every center, Gauss, tying, and
committed recovery location required by Section 9.2, the following quantities must
be finite and satisfy

$$
\boxed{
\|\mathbf A_\xi\times\mathbf A_\eta\|>0,
\qquad J>0.}
$$

The covariant vectors and reciprocal bases must also be finite. Duplicate nodes,
self-intersection, zero-area mappings, and nonpositive determinants fail before
stiffness or recovery is committed. No failed location is omitted, averaged, clamped,
or repaired. This feature defines no calibrated distortion, aspect, warpage, or
director-angle threshold; `NR-O03` and `NR-O04` are not acceptance tests.

## 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 Fixed drilling stabilization

Let `R` be the index set of the eight director-tangent rotational coordinates in the
physical local ordering `q_20`. Form the finite strictly positive diagonal inventory

$$
\mathcal R_+=\{(K_{20})_{ii}\mid i\in\mathcal R,
(K_{20})_{ii}>0,\ (K_{20})_{ii}\text{ finite}\}.
$$

Every member has rotational-stiffness units `force*length`. For an otherwise valid
element, an empty `R+` is a deterministic numerical-validation failure. Define

$$
\boxed{k_{ref}=\min\mathcal R_+,\qquad k_d=10^{-3}k_{ref}},
$$

and use the same positive scalar at all four local drilling coordinates:

$$
\boxed{\mathbf K_d^l=k_d\mathbf I_4}.
$$

Translational diagonals have units `force/length` and shall never enter `R+`.
Off-diagonal entries, nonpositive entries, and nonfinite entries also do not enter
the minimum. The fixed coefficient is a project numerical-stability choice informed
by the thesis rule; it is not a physical constitutive parameter, an Abaqus algorithm,
or a claim of coefficient optimality.

There is no `rho_d`, coefficient sweep, stable-plateau selection, conditioning
calibration, artificial-energy ratio, or drilling-specific output contract in this
feature. Verification checks only the exact selection rule, dimensions, symmetry,
positivity, four-mode regularization, deterministic repeatability, and separation
from physical recovery.

### 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.
The positive uniform `k_d` block should remove only those drilling modes,
giving expected stabilized rank 18 and nullity 6. These are verification targets,
not substitutes for the normalized rigid-action checks defined in Section 17.1.

### 12.4 Energy identity

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`. `E_drill` is the internal quadratic identity associated
with `K_drill`; it is not a physical shell energy and is not emitted as a required
result. No drilling-energy ratio or warning threshold is defined.

### 12.5 DOF scaling for rank and conditioning evidence

Raw shell stiffness matrices mix translational and rotational units and therefore
must not be used for singular-value, eigenvalue, rank, or condition-number
acceptance. Define diagnostic coordinate scalings

$$
\mathbf q_{20}=\mathbf S_{20}\widehat{\mathbf q}_{20},\qquad
\mathbf S_{20}=\operatorname{blockdiag}_{I=1}^{4}(L_e\mathbf I_3,\mathbf I_2),
$$

$$
\mathbf q_g=\mathbf S_{24}\widehat{\mathbf q}_g,\qquad
\mathbf S_{24}=\operatorname{blockdiag}_{I=1}^{4}(L_e\mathbf I_3,\mathbf I_3).
$$

The corresponding scaled element matrices are

$$
\boxed{\widehat{\mathbf K}_{20}=\mathbf S_{20}^T\mathbf K_{20}\mathbf S_{20},
\qquad \widehat{\mathbf K}_e=\mathbf S_{24}^T\mathbf K_e^{24}\mathbf S_{24}}.
$$

Every entry of these matrices has units `force*length`. Rigid and other diagnostic
vectors are transformed by `q_hat=S^{-1}q`. Scaling is used only for algebraic
evidence; it does not modify the physical assembly, prescribed values, or solve.

For a global shell model, use the deterministic model length

$$
L_m=\max_{e\in\mathcal E_{active}}L_e
$$

and construct the analogous full-model scaling with `(L_m I_3,I_3)` per six-DOF
node. Restrict it to the stable free-DOF order as `S_f` and define

$$
\widehat{\mathbf K}_{ff}=\mathbf S_f^T\mathbf K_{ff}\mathbf S_f.
$$

Optional global condition and numerical-rank diagnostics use `K_hat_ff`; a valid `0 x 0 Kff`
case is classified separately and is not reported as singular. Numerical
condition calibration is not an implementation gate, and no raw mixed-unit matrix
may be used for any reported spectrum or condition estimate.

## 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. The
> equations through Section 15.3 live only in a 20-coordinate physical director
> chart; they are not a complete global 24-DOF nonlinear element.

### 15.1 Total Lagrangian kinematics

Let `q_20` be the physical nonlinear coordinate vector with three translations and
two director parameters per node. In the reference configuration,

$$
\mathbf x(\xi,\eta,\zeta;\mathbf q_{20})
=\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 chart coordinates `q_20,A` and `q_20,B`, define

$$
\mathbf h_A=\frac{\partial\mathbf x}{\partial q_{20,A}},
\qquad
\mathbf h_{AB}=\frac{\partial^2\mathbf x}
{\partial q_{20,A}\partial q_{20,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_{20})=
\frac{\partial\overline{\mathbf e}}{\partial q_{20,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_{20})
=\int_{\Omega_0}\mathbf B_A^T\mathbf s\,dV}.
$$

Collecting these entries gives the physical chart residual
`r_20^phys(q_20) in R^20`. A total 24-DOF residual is deliberately not written at
this point: drilling and external load can be combined with the physical residual
only after the nonlinear 20-to-24 coordinate map in Section 15.4 is selected.

### 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_{20,A}\partial q_{20,B}}.
$$

The consistent `20 x 20` physical-chart 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. These equations define
`K_20^phys=K_20^mat+K_20^geo`; they do not yet define congruence to a global
24-DOF tangent. For a constant generalized dead-load vector in this selected
20-coordinate chart, its chart load tangent is zero. 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 Unresolved 20-to-24 mapping and drilling boundary

Let a future finite global six-DOF coordinate vector be `q_g` and suppose an
approved nonlinear map exists:

$$
\mathbf q_{20}=\boldsymbol\Phi(\mathbf q_g),\qquad
\mathbf A(\mathbf q_g)=\partial\boldsymbol\Phi/\partial\mathbf q_g.
$$

Here `A` is a `20 x 24` mapping Jacobian. Conditional on that map, virtual work gives

$$
\mathbf r_g^{phys}=\mathbf A^T\mathbf r_{20}^{phys}.
$$

Its consistent global physical tangent must contain both the congruence and the
coordinate-map curvature term:

$$
\boxed{\mathbf K_g^{phys}=\mathbf A^T\mathbf K_{20}^{phys}\mathbf A+\sum_{a=1}^{20}r_{20,a}^{phys}(\partial^2\Phi_a/\partial\mathbf q_g^2)}.
$$

The candidate two-parameter director chart in Section 15.1 does not by itself
define `Phi`: a finite three-component global rotation coordinate, its
director-parallel gauge, chart update/recentering, and the first and second
derivatives of the complete map are still unspecified.

Likewise, a future drilling potential would have to be objective and expressed in
the same global coordinates:

$$
\Pi_d=\frac12\boldsymbol\gamma(\mathbf q_g)^T
\mathbf K_d^l(\mathbf q_g)\boldsymbol\gamma(\mathbf q_g),
\qquad
\mathbf r_g^{drill}=\frac{\partial\Pi_d}{\partial\mathbf q_g},
\qquad
\mathbf K_g^{drill}=\frac{\partial^2\Pi_d}{\partial\mathbf q_g^2}.
$$

A constant initial-frame diagonal penalty is not automatically objective under
large rotation. The nonlinear frame, drilling coordinate, scale update, and all
their first/second derivatives remain unresolved. Until those choices and
configuration-dependent load work are approved, neither a complete global residual
nor a complete global 24-DOF tangent exists. This does not block the current linear
kernel but blocks geometrically nonlinear implementation.

### 15.5 Newton equation and linear limit

Only after Section 15.4 is closed may the global residual and tangent be defined as

$$
\mathbf r_g=\mathbf A^T\mathbf r_{20}^{phys}+\mathbf r_g^{drill}-\mathbf f_{ext},
\qquad \mathbf K_T=\mathbf K_g^{phys}+\mathbf K_g^{drill}-\mathbf K_g^{load}.
$$

The conditional Newton equation at iteration `k` is

$$
\mathbf K_T(\mathbf q_g^{(k)})\Delta\mathbf q_g=-\mathbf r_g(\mathbf q_g^{(k)}),
\qquad \mathbf q_g^{(k+1)}=\mathbf q_g^{(k)}+\Delta\mathbf q_g.
$$

At the undeformed, stress-free state, recovering the current linear formulation
requires all of the following identities:

$$
\mathbf A_0=\mathbf T_p,\qquad \mathbf r_{20}^{phys}=\mathbf0,
\qquad \mathbf K_{20}^{geo}=\mathbf0.
$$

Then the coordinate-map curvature term vanishes and

$$
\mathbf K_g^{phys}=\mathbf T_p^T\mathbf K_{20}\mathbf T_p=\mathbf K_{phys}^{24}.
$$

If the future map also recovers `gamma=T_d q_g` and the approved constant linear
drill block in this limit, the complete tangent reduces to `K_e^24`. This is a
conditional consistency requirement, not evidence that a nonlinear `Phi` or
objective drilling potential has already been selected.

## 16. Algorithm pseudocode

### 16.1 Deterministic preprocessing

```text
for each shell element in stable source order:
    validate four distinct source nodes and source-order geometry
    compute L_e from the four consecutive midsurface edges
    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 any nonpositive pairwise incident-normal dot product before averaging
    reject nonfinite or zero incident normals and averaged vectors
    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

```text
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 the complete center/Gauss/tying/recovery geometry inventory
validate finite nonzero surface area and finite positive J at every required location
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
        reuse the validated point geometry and reciprocal bases
        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
collect positive finite tangent-rotation diagonals R+ from K20
require R+ nonempty; k_d = 1e-3 * min(R+); Kd_local = k_d * I4
Kdrill24 = T_d^T * Kd_local * T_d
Ke24 = Kphys24 + Kdrill24
form S20, S24, Khat20, and Khat_e for normalized algebraic evidence only
fint24 = Ke24 * q_g
residual24 = fint24 - f_CLOAD
return matrices, residual, transforms, and frames
```

### 16.3 Global linear-static lifecycle

```text
assemble all Ke24 contributions with stable element-local COO ordering
form model-length DOF scaling and Khat_ff for optional normalized diagnostics only
partition full K into Kff, Kfc, Kcf, Kcc in stable free/constrained order
factorize Kff before load assembly
assemble and deterministically aggregate nodal CLOAD
accept an exact-zero nodal moment separately; otherwise require rho_M <= 1e-12
solve Kff * df = Ff - Kfc * dc
reconstruct full displacement d
compute full residual r = K*d - F
recover shell rows and physical shell energy in stable source order
validate complete finite candidate state/output, then commit
```

### 16.4 Recovery

```text
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 required physical shell energy; emit no drilling-specific result
```

### 16.5 Future nonlinear tangent check

```text
given a future approved 20-coordinate finite-director chart:
    evaluate physical residual r20(q20) and consistent K20(q20)
    verify its directional derivative
given an approved global Phi map, objective drill potential, and load work:
    evaluate global residual rg(qg) and complete consistent K_T(qg)
    for several normalized perturbation directions p and decreasing h:
        compare K_T*p with [rg(qg+h*p)-rg(qg-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 nonzero surface area and finite positive `J` at every center, Gauss, tying,
  and committed recovery point.
- Scaled `K20`, `Kphys24`, `Kdrill24`, and `Ke24` symmetry and spectrum.
- Transformation work/energy invariance.
- Physical rigid modes satisfy normalized scaled stiffness action.
- Stabilized free-element nullity is exactly six; accepted non-rigid physical modes
  have positive physical energy.
- Pure drill vectors have zero physical-shell energy and positive action under the
  fixed numerical drilling block.
- Consistent force/length unit rescaling leaves dimensionless decisions unchanged.

For any nonzero scaled stiffness under test, the approved normalized checks are
the following. Here `r_hat` denotes a constructed scaled rigid-mode vector, not the
assembled residual.

$$
e_{sym}=\|\widehat{\mathbf K}-\widehat{\mathbf K}^T\|_F/\|\widehat{\mathbf K}\|_F\le10^{-12},
$$

$$
e_{rigid}=\|\widehat{\mathbf K}\widehat{\mathbf r}\|_2/(\|\widehat{\mathbf K}\|_2\|\widehat{\mathbf r}\|_2)\le10^{-10},
$$

$$
e_{frame}=\|\mathbf R^T\mathbf R-\mathbf I\|_F\le10^{-12}.
$$

For a nonzero transformation-energy case,

$$
e_T=|E_g-E_l|/(|E_g|+|E_l|)\le10^{-12}.
$$

A zero matrix/vector denominator is invalid test construction and is not clamped to
pass. When both transformation energies are exactly zero, the rigid or pure-null
case is classified by its separate stiffness-action test. Linear-system residual
and global-equilibrium evidence retain the normalized `1e-10` requirement.

### 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 Nonblocking accuracy studies

Thin/thick, distorted, warped, pinched-cylinder, NAFEMS LE3 and Scordelis-Lo studies
remain useful for documenting the original MITC4 element's known shear- and
membrane-locking limits. They are not required implementation-completion tests and
do not authorize MITC4+ or an expanded input subset. `NR-O03` and `NR-O04` are
explicitly removed from the acceptance scope.

### 17.4 Fixed drilling checks

For representative valid element kernels, verify:

- `R+` contains only finite positive physical tangent-rotation diagonals;
- `k_d=10^-3 min(R+)` and `K_d^l=k_d I4` exactly;
- `K_drill^24` is symmetric and positive on each pure drilling coordinate;
- four nonphysical drilling null modes are removed while the six physical rigid
  modes satisfy the normalized action test;
- physical generalized strain, resultant and stress recovery is unchanged by a
  pure drilling vector.

No coefficient sweep, plateau, condition threshold, response-sensitivity criterion,
or drilling-energy warning is part of this check.

### 17.5 Reference-comparison boundary

Abaqus comparisons block only on matched global `U1/U2/U3` rows under
`tolerance_c=1e-9+1e-6*reference_scale_c`, where `reference_scale_c` is the maximum
absolute finite Abaqus value for the same component. `UR1/UR2/UR3` uses the same
formula but an exceedance emits only a deterministic nonblocking warning.
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, the first/second derivatives of `Phi`, the coordinate-map curvature
term, objective drilling, 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 patch/reference evidence; broader convergence is a known limitation study |
| Membrane locking on distorted curved meshes | slow or nonuniform convergence | document original MITC4 limitation; optional later studies do not alter the current gate |
| 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 | finite nonzero surface area and finite positive `J` at every required location |
| 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 |
| Invalid drilling reference inventory | absent or nonpositive numerical regularization | require nonempty finite positive physical rotational diagonals and fail deterministically otherwise |
| Mixed-unit drilling scale | unit-dependent stabilization | exclude every translational diagonal; use only the physical tangent-rotation block |
| 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 |
| Missing nonlinear 20-to-24 map Hessian | incomplete global tangent | require `Phi`, `A`, and coordinate-map curvature term |
| 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; fixed positive rotational-diagonal stabilization | Numerical Review for formula consistency |
| `006-010` | isotropic plane stress, one centered constant-thickness layer | I/O validation |
| `011-016` | pairwise orientation, deterministic averaging/frames, basic finite/positive geometry predicates | I/O validation |
| `017-020` | global 24-DOF virtual work; `rho_M <= 1e-12`; distributed loads excluded | I/O diagnostic/schema |
| `024-029` | deterministic element buffers, partitioned linear lifecycle, full-residual reaction | planning |
| `039-048` | nodal/global and shell/local recovery inventory, units, identities, physical shell energy | I/O schema |
| `049-057` | normalized invariants, patches, fixed drilling, declared references and equilibrium | Numerical Review/reference/physics |
| `058-064` | exact B33 mixed tolerance; U blocking and UR warning-only | reference verification |
| `065-072` | exact existing S4/S4R paths, immutability and displacement-only boundary | reference model |

### 20.1 Numerical Review revision traceability

The first review findings map to this revision as follows.

`NR-C01` maps to Sections 9.2-9.3, 16.2, and 17.1/17.3. Basic point inventory and
finite/positive predicates remain; `NR-O03`/`NR-O04` calibration is removed.

`NR-C02` maps to Sections 12.2 and 17.4. The project decision replaces candidate
normalization and plateau work with the exact fixed rotational-diagonal rule.

`NR-C03` maps to Sections 12.3/12.5 and 17.1. Element/global DOF scaling is fixed
while the numerical-rank threshold remains open.

`NR-C04` maps to Sections 7.1-7.2, which fix the global 24-DOF weak form and
transpose notation. `NR-C05` maps to Sections 15.1-15.5 and 17.6, which expose the
physical-chart/global-map boundary and required map Hessian term.

`NR-D01` maps to Sections 6.2 and 16.3, which fix the exact-zero case and
`rho_M <= 1e-12`. `NR-D02` maps to Sections 12.5 and 17.1, which fix the approved
normalized algebraic metrics.

## 21. Open issues and downstream handoff

### 21.1 Current numerical-review boundary

No calibration decision remains open for the linear implementation. The first review
approved `rho_M <= 1e-12`, the normalized algebraic checks, MITC tying/component signs,
constitutive law, quadrature, and recovery signs. The approved policy fixes drilling
and U/UR tolerance and removes drilling-energy calibration plus `NR-O03`/`NR-O04`.
Numerical Review shall now check internal consistency and may not treat those removed
items or an expanded reference portfolio as missing evidence.

### 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 orientation, basic topology/Jacobians, unsupported loads, and recovery;
  encode exact-zero nodal moment separately and enforce `rho_M <= 1e-12` without
  a denominator clamp.
- 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 physical shell energy. Do not add drilling-specific datasets.
- Preserve local frame and natural-coordinate identity without location averaging.

### 21.3 Reference Model handoff

- Record the existing `reference/shell/` S4 and `reference/shellR/` S4R input and
  displacement CSV paths, but compare only formulation-independent global displacement
  evidence rather than claiming element equivalence.
- Use the exact B33 mixed tolerance; do not add administrative metadata or portfolio gates.
- Do not create, repair, rename, or run reference artifacts during this formulation
  gate.

### 21.4 Implementation Planning handoff

- Begin after the Numerical Review, I/O contract, and lightweight Reference Model
  inventory agree with this formulation.
- Translate the deterministic preprocessing, 24-to-20 transform, tying projection,
  fixed physical/drill split, quadrature, recovery, and required invariants into
  Harness Step drafts with `RED -> GREEN -> VERIFY` tests before production changes.
- Obtain user approval of the multi-Step draft before creating phase-planning files;
  Harness execution requires a separate explicit user request.
- 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, the complete `Phi: R24 -> R20` map and its first/second
derivatives, 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 physical-chart residual/tangent and the conditional global
mapping identity but deliberately does not close those product decisions.