FESADev/docs/nonlinear-static-newton-raphson/formulation.md

# Nonlinear Static Newton–Raphson 정식화

## 1. Metadata

- feature_id: nonlinear-static-newton-raphson
- source_requirement: docs/nonlinear-static-newton-raphson/requirements.md
- source_requirement_status: approved
- source_requirement_approval_date: 2026-08-23
- source_requirement_sha256: 910E351F3F90FFA42F07D15206B55E62C05EAE821E250424430148337506EBF2
- source_research: docs/nonlinear-static-newton-raphson/research.md
- source_research_status: passed
- source_research_sha256: C920A6B02381230D1DB419D73FE767D611B36720D27D6410507BA88C3CA28EEE
- status: ready-for-numerical-review
- owner_agent: formulation-agent
- date: 2026-08-23
- formulation_scope: fixed-initial-reference Total Lagrangian, materially linear and
  geometrically nonlinear original-style FESA-MITC4, Full Newton
- independence: FESA는 독립 solver이며 Abaqus, Midas 또는 다른 solver의 요소, 회전,
  stabilization, 적분 또는 상태 갱신 동등성을 주장하지 않는다.

이 문서는 Stage 3의 수학·알고리즘 계약이다. C++ API, parser, HDF5 path/schema, exact
reference artifact, benchmark blocking 분류 및 수치 검토 tolerance는 정의하지 않는다.
기존 docs/linear-static-mitc4-shell/formulation.md의 §1–14에서 승인된 node order,
initial director, local sign, MITC tying, constitutive boundary, 2 x 2 x 2 integration과
physical recovery를 보존한다. 그 문서 §15는 gap inventory와 conditional identity로만
사용하며 이 문서의 nonlinear 식을 대신하지 않는다.

## 2. Scope and Assumptions

### 2.1 Analysis boundary

- 해석 종류는 정확히 하나의 nonlinear static step이다.
- 하중 제어는 dimensionless scalar factor \(\lambda\in[0,1]\)의 fixed increment path다.
- 각 증분은 Full Newton으로 평형을 찾는다. Modified Newton, tangent reuse, line search,
  damping, cutback, retry 및 arc-length는 없다.
- 모든 residual, tangent, load-work derivative와 convergence quantity는 같은 current
  trial identity에서 평가한다.
- formulation은 한 고정 initial reference configuration \(\Omega_0\)만 사용한다.
  accepted increment 뒤에도 reference geometry를 재설정하지 않는다.
- inertia, mass, damping 및 time integration은 N/A다.

### 2.2 Element and material boundary

- Source S4와 S4R은 같은 four-node original-style FESA-MITC4를 선택한다.
- Element는 Reissner–Mindlin-type degenerated continuum shell이며, thickness stretch와
  \(E_{33}\), \(S_{33}\)은 제외한다.
- Material은 centered homogeneous isotropic linear elasticity다.
- Geometry는 large displacement와 finite rotation을 허용하지만 material applicability는
  small elastic strain이다. Material history/state variable은 없다.
- MITC4+ membrane field, reduced integration, hourglass control, EAS, physical drilling
  constitutive law, follower load 및 distributed load는 포함하지 않는다.

### 2.3 Coordinates, signs and units

- Natural coordinates are
  \((\xi,\eta,\zeta)\in[-1,1]^3\).
- Node layout is
  \(1=(-1,-1)\), \(2=(+1,-1)\), \(3=(+1,+1)\), \(4=(-1,+1)\).
- \(+\zeta\) and the initial director define the positive shell face.
- Global nodal DOF order is
  \([UX,UY,UZ,URX,URY,URZ]\).
- Translation and nodal force are global and have dimensions length and force.
- Rotation coordinates/corrections are dimensionless. Nodal moments have dimension
  force-length.
- Constitutive strain, second Piola–Kirchhoff stress and physical shell recovery use
  the fixed initial local frame at each reference midsurface location. No rotating
  output frame is inserted into the material tangent.
- User input uses one consistent unit system. FESA does not infer or convert units.

## 3. Primary Variables, DOFs and State

### 3.1 Fixed reference data

At node \(I\), let

\[
\mathbf X_I\in\mathbb R^3,\qquad
\mathbf Q_I^0=
\begin{bmatrix}\mathbf a_I^0&\mathbf b_I^0&\mathbf d_I^0\end{bmatrix}
\in SO(3)
\]

be the initial midsurface position and the approved initial right-handed triad.
Initial directors and tangent frames are constructed exactly by the approved linear
MITC4 deterministic incident-element rules. Thickness \(t>0\), \(E>0\), and
\(-1<\nu<0.5\) are fixed.

### 3.2 Global nonlinear coordinates

The nodal nonlinear coordinate is

\[
\mathbf q_I=
\begin{bmatrix}\mathbf u_I\\\boldsymbol\phi_I\end{bmatrix}
\in\mathbb R^6.
\]

\(\mathbf u_I\) is global midsurface displacement. \(\boldsymbol\phi_I\) is a
fixed-initial, fixed-global-basis rotation vector, and

\[
\mathbf R_I(\boldsymbol\phi_I)
=\exp([\boldsymbol\phi_I]_\times),\qquad
\mathbf Q_I=\mathbf R_I\mathbf Q_I^0,\qquad
\mathbf d_I=\mathbf Q_I\mathbf e_3 .
\]

The mathematical branch is the principal exponential chart
\(\|\boldsymbol\phi_I\|<\pi\). A target or trial exactly on the non-unique
\(\pi\)-rotation branch, an ambiguous logarithm, or nonfinite rotation data fails
closed. Numerical branch-distance tolerance is owned by Numerical Review; this
formulation introduces no hidden clamp.

The 24-vector for one element is

\[
\mathbf q_g=
\begin{bmatrix}
\mathbf u_1&\boldsymbol\phi_1&
\cdots&
\mathbf u_4&\boldsymbol\phi_4
\end{bmatrix}^{T}\in\mathbb R^{24}.
\]

Constrained and free identities are the fixed global components of this vector.

### 3.3 Multiplicative update and retraction

For an unchanged Newton correction
\(\Delta\boldsymbol\phi_I\), define

\[
\mathbf Q_I^+
=\exp([\boldsymbol\phi_I+\Delta\boldsymbol\phi_I]_\times)\mathbf Q_I^0
=\Delta\mathbf R_I\,\mathbf Q_I,
\]

\[
\boxed{
\Delta\mathbf R_I=
\exp([\boldsymbol\phi_I+\Delta\boldsymbol\phi_I]_\times)
\exp(-[\boldsymbol\phi_I]_\times)\in SO(3).}
\]

Thus the update is left multiplication in the spatial/global frame. The correction
is dimensionless, expressed in the fixed global basis, and is applied without a
scalar modification. The authoritative orientation is \(\mathbf Q_I\); after every
update \(\boldsymbol\phi_I\), \(\mathbf Q_I\), and \(\mathbf d_I\) must satisfy the
boxed identity. Corrections are not multiplied successively as independent
incremental rotation tensors; the fixed-reference exponential coordinate is updated
first and \(\Delta\mathbf R_I\) is then derived. Therefore a final
\(\boldsymbol\phi_I\) defines one path-independent triad inside the principal chart.

Rodrigues' map is

\[
\exp([\boldsymbol\phi]_\times)
=\mathbf I+\frac{\sin\theta}{\theta}[\boldsymbol\phi]_\times
+\frac{1-\cos\theta}{\theta^2}[\boldsymbol\phi]_\times^2,
\qquad \theta=\|\boldsymbol\phi\|,
\]

with its analytic series at \(\theta=0\). No normalization of \(\mathbf Q\) is part
of the mathematical update. A finite implementation may apply a checked
orthogonal projection only as a representation safeguard if Numerical Review
approves its invariant tolerance; it must not change residual/tangent state identity.

### 3.4 Exact first and second rotation variations

For a direction \(\mathbf p\), define the left Jacobian

\[
\mathbf J_l(\boldsymbol\phi)=
\mathbf I+
\frac{1-\cos\theta}{\theta^2}[\boldsymbol\phi]_\times+
\frac{\theta-\sin\theta}{\theta^3}[\boldsymbol\phi]_\times^2 .
\]

Then the spatial infinitesimal rotation induced by the additive global-coordinate
variation \(\mathbf p\) is

\[
\boldsymbol\omega_{\mathbf p}
=\mathbf J_l(\boldsymbol\phi)\mathbf p,
\qquad
D\mathbf d[\mathbf p]
=\boldsymbol\omega_{\mathbf p}\times\mathbf d .
\]

For directions \(\mathbf p,\mathbf s\),

\[
\boxed{
D^2\mathbf d[\mathbf p,\mathbf s]
=D\mathbf J_l(\boldsymbol\phi)[\mathbf s]\mathbf p\times\mathbf d
+\boldsymbol\omega_{\mathbf p}\times
(\boldsymbol\omega_{\mathbf s}\times\mathbf d).}
\]

The symmetric value is obtained by the analytically identical exchange
\(\mathbf p\leftrightarrow\mathbf s\); disagreement is an implementation error.
\(D\mathbf J_l\) is the direct derivative of the displayed analytic expression,
including scalar coefficient derivatives and the series branch. It is not omitted.

Equivalently, all required derivatives may be evaluated through the exact Fréchet
primitives

\[
D\exp_{\mathbf A}[\mathbf H]
=\int_0^1
\exp((1-s)\mathbf A)\mathbf H\exp(s\mathbf A)\,ds,
\]

\[
D^2\exp_{\mathbf A}[\mathbf H,\mathbf K]
=D\!\left(D\exp_{\mathbf A}[\mathbf H]\right)[\mathbf K],
\]

where the second expression is evaluated symmetrically. These definitions, the
Rodrigues form, and the directional-derivative verification must give the same
first/second variation.

### 3.5 Trial and committed state

A complete state identity contains at least

\[
\mathcal S=
\{\lambda,\mathbf u,\boldsymbol\phi,\mathbf Q,\mathbf d,
\mathbf f_{int},\mathbf f_{ext},\mathbf r,
\text{reaction},\text{physical recovery rows}\}.
\]

Committed and trial instances never alias. Element drilling quantities are derived
from the same trial and introduce no material history. A successful increment
atomically transfers the entire validated trial identity. Nonconvergence, nonfinite
data, invalid geometry/rotation, solve failure, recovery failure, or output failure
discards the candidate and preserves the preceding committed state.

## 4. Strong Form and Boundary Conditions

### 4.1 Total Lagrangian parent continuum

On the fixed reference volume \(\Omega_0\),

\[
\operatorname{Div}_0\mathbf P+\mathbf B_0=\mathbf0,\qquad
\mathbf P=\mathbf F\mathbf S ,
\]

\[
\mathbf F=\frac{\partial\mathbf x}{\partial\mathbf X},\qquad
\mathbf E=\frac12(\mathbf F^T\mathbf F-\mathbf I),\qquad
\mathbf S=\mathbb C_{ps}:\overline{\mathbf E}.
\]

The parent boundary is split into \(\Gamma_u\cup\Gamma_t\):

\[
\mathbf x=\overline{\mathbf x}\quad\text{on }\Gamma_u,\qquad
\mathbf P\mathbf N_0=\overline{\mathbf T}_0\quad\text{on }\Gamma_t .
\]

For this executable subset \(\mathbf B_0=\mathbf0\) and distributed
\(\overline{\mathbf T}_0=\mathbf0\). External work is supplied only by approved
nodal CLOAD and nodal prescribed global coordinates. The continuum natural term
documents the lineage but does not authorize pressure, edge traction, body force,
or follower loading.

### 4.2 Translational prescribed path

For a constrained translation component \(c\), step-start value
\(u_c^s\), factor-one absolute target \(u_c^\star\), and current factor \(\lambda\),

\[
\boxed{u_c(\lambda)=u_c^s+\lambda(u_c^\star-u_c^s).}
\]

The step has no prior history in this release, but the formula preserves the
step-start identity explicitly.

### 4.3 Finite rotational prescribed path

For a constrained global rotation-vector component \(c\),

\[
\boxed{
\phi_c(\lambda)=\phi_c^s+
\lambda(\phi_c^\star-\phi_c^s).}
\]

\(\boldsymbol\phi^\star\) is the factor-one absolute global rotation-vector target,
not a Newton correction. At an increment target, constrained coordinates are set
to the boxed absolute values before the first residual evaluation. Every Newton
solve then uses

\[
\boxed{\Delta\phi_c=0.}
\]

Thus target advance and Newton correction are distinct. The triad path is
\(\mathbf Q(\lambda)=\exp([\boldsymbol\phi(\lambda)]_\times)\mathbf Q^0\) only for
components fixed by boundary data; free components remain solution variables.
This componentwise path is deterministic, reaches the exact source target at
\(\lambda=1\), and introduces no follower or controlled-displacement semantics.

### 4.4 Nodal load path and drilling-moment rejection

Approved nodal force and immutable source-moment coefficients
\(\mathbf F_I^\star,\mathbf M_{I,src}^\star\) remain fixed in global axes:

\[
\mathbf F_I(\lambda)=\lambda\mathbf F_I^\star,\qquad
\mathbf M_{I,src}(\lambda)=\lambda\mathbf M_{I,src}^\star .
\]

The factor changes magnitude only. A zero aggregated nodal moment is accepted as a
separate exact-zero case. Every nonzero moment is checked against the same current
trial director used by the physical residual:

\[
\boxed{
\rho_{M,I}=
\frac{|\mathbf d_I\cdot\mathbf M_{I,src}^\star|}
{\|\mathbf M_{I,src}^\star\|}
\le 10^{-12}.}
\]

Failure at preprocessing or at any entered trial is
unsupported-drilling-load and occurs before that trial residual/tangent can be
accepted. The denominator is never clamped. This admissibility check does not
project, rotate, shorten, overwrite, or otherwise mutate the source vector; source
identity and the proportional fixed-global path remain exact.

The shell has no physical director-parallel work channel. After admissibility, the
exact covector representing the immutable source moment on the current physical
two-DOF virtual-rotation space is

\[
\boxed{
\mathbf M_{I,phys}(\mathbf q;\lambda)
=\lambda\mathbf P_I(\mathbf q)\mathbf M_{I,src}^\star,\qquad
\mathbf P_I=\mathbf I-\mathbf d_I\mathbf d_I^T.}
\]

Thus an admitted source component within the approved \(10^{-12}\) direction band
is removed from physical work exactly and is not sent to the numerical drill
channel. \(\mathbf M_{I,phys}\) is an effective work covector, not a mutated source
load. Its state dependence is solely the restriction to the changing admissible
director-tangent space; it does not rotate the fixed-global source direction and
does not introduce follower-load semantics. Residual and tangent must include
\(D\mathbf P_I\) exactly as specified in §11.2.

## 5. Weak or Variational Form

Let \(\mathcal V_{24}(\lambda)\) be the admissible trial manifold satisfying the
absolute constraints and let \(\mathcal V_{24}^0\) contain variations with zero
constrained global components. The equilibrium virtual-work statement is

\[
\boxed{
\delta\mathcal G(\mathbf q_g;\lambda)
=\delta\Pi_{phys}
+\delta\Pi_{drill}
-\delta W_F
-\delta W_M=0
\quad\forall\delta\mathbf q_g\in\mathcal V_{24}^0 .}
\]

The physical internal virtual work is

\[
\delta\Pi_{phys}
=\int_{\Omega_0}
\delta\overline{\mathbf e}^{T}\mathbf s\,dV_0,
\qquad
\mathbf s=\mathbf C_5\overline{\mathbf e}.
\]

The translational dead-load work is

\[
W_F=\lambda\sum_I\mathbf F_I^\star\cdot\mathbf u_I .
\]

For a virtual rotation-coordinate variation
\(\delta\boldsymbol\phi_I\), define the induced spatial rotation and its physical
director-tangent part by

\[
\delta\boldsymbol\vartheta_I
=\mathbf J_l(\boldsymbol\phi_I)\delta\boldsymbol\phi_I,\qquad
\delta\boldsymbol\vartheta_I^{phys}
=(\mathbf I-\mathbf d_I\mathbf d_I^T)
\delta\boldsymbol\vartheta_I
=\mathbf d_I\times\delta\mathbf d_I .
\]

\[
\boxed{
\delta W_M
=\sum_I\mathbf M_{I,phys}(\mathbf q;\lambda)\cdot
\delta\boldsymbol\vartheta_I
=\lambda\sum_I\mathbf M_{I,src}^\star\cdot
\delta\boldsymbol\vartheta_I^{phys}.}
\]

Here \(\mathbf M_{I,src}^\star\) is the immutable global source vector and
\(\mathbf M_{I,phys}=\lambda\mathbf P_I\mathbf M_{I,src}^\star\) is the exact
effective physical-work covector from §4.4. The two expressions are identical
because \(\delta\boldsymbol\vartheta_I^{phys}=\mathbf P_I
\delta\boldsymbol\vartheta_I\). Admissibility permits a source component within
the approved ratio band; it does not assert
\(\mathbf M_{I,src}^\star\cdot\mathbf d_I=0\) exactly. Projection removes that
component from physical work, and the numerical director-parallel rotation does no
external work. Configuration dependence of the generalized load arises from the
exact rotation-coordinate and projector derivatives; the immutable source vector
itself is never projected in storage, rotated, or shortened.

The physical shell and drilling terms are conservative inside their valid charts.
An immutable constant-spatial source moment restricted to the physical tangent
space as above defines a work one-form on finite rotations, not a global scalar
rotation potential. Its exact load-work tangent is therefore generally nonsymmetric
even though it contains no follower direction update. With no applied moment, the
complete unreduced internal/drill tangent is a symmetric Hessian in exact arithmetic.
Numerical symmetry expectations and norms belong to Numerical Review.

## 6. Discretization

### 6.1 Bilinear interpolation

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

\[
\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}
\]

Mandatory algebraic checks are

\[
\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.
\]

These checks apply independently of deformation.

### 6.2 Reference and current shell geometry

With \(z=t\zeta/2\),

\[
\boxed{
\mathbf X(\xi,\eta,\zeta)
=\sum_I N_I\mathbf X_I+
z\sum_I N_I\mathbf d_I^0,}
\]

\[
\boxed{
\mathbf x(\xi,\eta,\zeta;\mathbf q_g)
=\sum_I N_I(\mathbf X_I+\mathbf u_I)+
z\sum_I N_I\mathbf d_I.}
\]

Straight material lines through thickness remain straight. The current director is
unit length by the rotation map. Thickness is constant and is not a state variable.

The continuum director field used by geometry is

\[
\mathbf D(\xi,\eta)=\sum_I N_I\mathbf d_I .
\]

It is bilinearly interpolated and is not normalized before entering
\(\mathbf x\), \(\mathbf g_\alpha\), strain, tying, residual, or tangent. This is
the selected original degenerated-continuum convention; silently replacing it by
\(\mathbf D/\|\mathbf D\|\) changes the kinematics. Its exact variations are simply

\[
D\mathbf D=\sum_I N_I D\mathbf d_I,\qquad
D^2\mathbf D=\sum_I N_I D^2\mathbf d_I.
\]

When a unit interpolated director is needed only for a checked diagnostic/current
frame, use \(\overline{\mathbf d}=\mathbf D/\|\mathbf D\|\) after requiring
\(\|\mathbf D\|>0\). With
\(\mathbf n=\mathbf v/\|\mathbf v\|\) and
\(\mathbf P=\mathbf I-\mathbf n\mathbf n^T\),

\[
D\mathbf n[\mathbf p]=\frac{\mathbf P\mathbf p}{\|\mathbf v\|},
\]

\[
\boxed{
D^2\mathbf n[\mathbf p,\mathbf s]
=-\frac{
(\mathbf n\cdot\mathbf s)\mathbf P\mathbf p+
(\mathbf n\cdot\mathbf p)\mathbf P\mathbf s+
\mathbf n(\mathbf p\cdot\mathbf P\mathbf s)}
{\|\mathbf v\|^2}.}
\]

No normalized diagnostic frame feeds back into the physical constitutive law.
At every trial require

\[
\mathbf d_I\cdot\mathbf d_J>-1\quad(I<J),
\qquad
\|\mathbf D(\xi,\eta)\|>0
\]

at every required tying, Gauss, center and recovery location. Thus an exactly
opposed nodal pair or a degenerate interpolated director fails closed rather than
choosing a replacement axis. Scale-aware proximity tolerances, if required, belong
to Numerical Review.

Define

\[
\mathbf G_\alpha=\mathbf X_{,\alpha},\qquad
\mathbf g_\alpha=\mathbf x_{,\alpha},\qquad
\alpha\in\{\xi,\eta,\zeta\}.
\]

Constant translation is reproduced exactly. For any proper rigid
\(\mathbf R\) and translation \(\mathbf c\),

\[
\mathbf x_I=\mathbf R\mathbf X_I+\mathbf c,\qquad
\mathbf Q_I=\mathbf R\mathbf Q_I^0
\]

gives \(\mathbf g_\alpha=\mathbf R\mathbf G_\alpha\) and zero physical strain.

## 7. Kinematics and the Nonlinear 24-to-20 Map

### 7.1 Fixed-reference strain

\[
\mathbf F=\mathbf g_\alpha\otimes\mathbf G^\alpha,\qquad
\mathbf G^\alpha\cdot\mathbf G_\beta=\delta^\alpha_\beta,
\]

\[
\boxed{
E_{\alpha\beta}
=\frac12(\mathbf g_\alpha\cdot\mathbf g_\beta
-\mathbf G_\alpha\cdot\mathbf G_\beta).}
\]

This is Green–Lagrange strain. Second Piola–Kirchhoff stress is its
work-conjugate stress.

### 7.2 Exact physical director chart

Let

\[
\mathbf c_I=\mathbf d_I^0\times\mathbf d_I,\qquad
s_I=\|\mathbf c_I\|,\qquad
c_I=\mathbf d_I^0\cdot\mathbf d_I.
\]

For \(c_I>-1\), define the minimal swing vector

\[
\boxed{
\boldsymbol\chi_I=
\begin{cases}
\mathbf0,&s_I=0,\ c_I=1,\\[1mm]
\dfrac{\operatorname{atan2}(s_I,c_I)}{s_I}\mathbf c_I,
&s_I>0 .
\end{cases}}
\]

\(\boldsymbol\chi_I\cdot\mathbf d_I^0=0\) and
\(\exp([\boldsymbol\chi_I]_\times)\mathbf d_I^0=\mathbf d_I\).
The antipodal state \(c_I=-1\) is not a unique two-DOF director chart and fails
closed.

Define the two physical rotation coordinates

\[
\boldsymbol\chi_I^{2}=
\begin{bmatrix}
\mathbf a_I^0\cdot\boldsymbol\chi_I\\
\mathbf b_I^0\cdot\boldsymbol\chi_I
\end{bmatrix}\in\mathbb R^2
\]

and

\[
\boxed{
\mathbf q_{20}=\boldsymbol\Phi(\mathbf q_g)
=
\begin{bmatrix}
\mathbf u_1&\boldsymbol\chi_1^2&
\cdots&
\mathbf u_4&\boldsymbol\chi_4^2
\end{bmatrix}^{T}\in\mathbb R^{20}.}
\]

This map discards only the rotation about the director gauge. It is an exact finite
map, not the unresolved linear projection from the earlier gap inventory.

### 7.3 Map Jacobian and Hessian

Define

\[
\mathbf A(\mathbf q_g)
=D\boldsymbol\Phi(\mathbf q_g)
\in\mathbb R^{20\times24},
\qquad
\mathbf H_a^\Phi(\mathbf q_g)
=D^2\Phi_a(\mathbf q_g)
\in\mathbb R^{24\times24}.
\]

The translation blocks are identity and have zero Hessian. Rotation blocks are
obtained by the exact chain

\[
\boldsymbol\phi
\xrightarrow{\exp}
\mathbf d
\xrightarrow{\text{minimal swing}}
\boldsymbol\chi
\xrightarrow{(\mathbf a^0,\mathbf b^0)}
\boldsymbol\chi^2 .
\]

For implementation, first use §3.4 for
\(D\mathbf d,D^2\mathbf d\), then differentiate
\(\mathbf c=\mathbf d^0\times\mathbf d\),
\(s=\|\mathbf c\|\), \(c=\mathbf d^0\cdot\mathbf d\), and
\(\alpha=\operatorname{atan2}(s,c)\). The exact primitives are

\[
D\|\mathbf v\|[\mathbf p]
=\frac{\mathbf v\cdot\mathbf p}{\|\mathbf v\|},
\]

\[
D^2\|\mathbf v\|[\mathbf p,\mathbf s]
=\frac{\mathbf p\cdot\mathbf s}{\|\mathbf v\|}
-\frac{(\mathbf v\cdot\mathbf p)(\mathbf v\cdot\mathbf s)}
{\|\mathbf v\|^3},
\]

\[
D\operatorname{atan2}(y,x)[p]
=\frac{x\,Dy[p]-y\,Dx[p]}{x^2+y^2},
\]

with the direct derivative of this last expression for the mixed second variation.
The analytic \(s\rightarrow0,c\rightarrow1\) series must be used; neither first nor
second derivative may be obtained by normalizing a zero cross product.

At the undeformed state,

\[
\boxed{
\mathbf A_0=\mathbf T_p,}
\]

where \(\mathbf T_p\) is the approved linear 20-by-24 physical transform, and

\[
D\boldsymbol\chi_I^2[\delta\boldsymbol\phi_I]
=
\begin{bmatrix}
(\mathbf a_I^0)^T\\
(\mathbf b_I^0)^T
\end{bmatrix}
\delta\boldsymbol\phi_I .
\]

The matrices \(\mathbf H_a^\Phi\) are generally nonzero away from the linear state.
They are required in the global tangent.

### 7.4 First and second strain variations

For a physical coordinate \(q_{20,A}\), 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}}.
\]

Then

\[
\boxed{
E_{\alpha\beta,A}
=\frac12(
\mathbf h_{A,\alpha}\cdot\mathbf g_\beta+
\mathbf g_\alpha\cdot\mathbf h_{A,\beta}),}
\]

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

Director-coordinate blocks retain the exact first and second derivative of
\(\exp([\boldsymbol\chi]_\times)\mathbf d^0\). Translation-translation
\(\mathbf h_{AB}\) is zero, but the quadratic gradient products remain and generate
stress-dependent geometric stiffness.

## 8. MITC4 Tying and Local Engineering Measures

### 8.1 Tying locations

\[
T_{\xi-}=(0,-1,0),\quad
T_{\xi+}=(0,+1,0),\quad
T_{\eta-}=(-1,0,0),\quad
T_{\eta+}=(+1,0,0).
\]

The first pair supplies \(E_{\xi\zeta}\) and the second pair supplies
\(E_{\eta\zeta}\).

### 8.2 Nonlinear assumed shear

\[
\boxed{
\overline E_{\xi\zeta}(\xi,\eta)
=\frac{1-\eta}{2}E_{\xi\zeta}^{DI}(T_{\xi-})
+\frac{1+\eta}{2}E_{\xi\zeta}^{DI}(T_{\xi+}),}
\]

\[
\boxed{
\overline E_{\eta\zeta}(\xi,\eta)
=\frac{1-\xi}{2}E_{\eta\zeta}^{DI}(T_{\eta-})
+\frac{1+\xi}{2}E_{\eta\zeta}^{DI}(T_{\eta+}).}
\]

All in-plane components remain direct. The same linear projection operator is
applied to the strain and both variations:

\[
\overline{\mathbf E}=\mathcal P_{MITC}\mathbf E,\qquad
D\overline{\mathbf E}=\mathcal P_{MITC}D\mathbf E,\qquad
D^2\overline{\mathbf E}=\mathcal P_{MITC}D^2\mathbf E .
\]

Tying geometry/directors and the trial identity must be identical for residual and
tangent. Direct shear in any one of the three is prohibited.

### 8.3 Fixed reference local frame

At each reference midsurface location, construct the approved initial frame
\((\mathbf e_1^0,\mathbf e_2^0,\mathbf e_3^0)\). Transform the assumed tensor into
that frame and define

\[
\boxed{
\overline{\mathbf e}=
\begin{bmatrix}
E_{11}&E_{22}&2E_{12}&2E_{13}&2E_{23}
\end{bmatrix}^{T}.}
\]

Let

\[
\mathbf B_A=\frac{\partial\overline{\mathbf e}}
{\partial q_{20,A}}\in\mathbb R^5,\qquad
\mathbf G_{AB}=
\frac{\partial^2\overline{\mathbf e}}
{\partial q_{20,A}\partial q_{20,B}}\in\mathbb R^5 .
\]

The fixed frame has no trial derivative. A rotating recovery frame, if ever
defined downstream, cannot be substituted into this constitutive calculation.

## 9. Constitutive Contract

Let

\[
G=\frac{E}{2(1+\nu)},\qquad
\kappa_s=\frac56,
\]

\[
\mathbf C_{ps}=
\frac{E}{1-\nu^2}
\begin{bmatrix}
1&\nu&0\\
\nu&1&0\\
0&0&(1-\nu)/2
\end{bmatrix},
\]

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

The second Piola–Kirchhoff engineering stress vector is

\[
\boxed{
\mathbf s=
\begin{bmatrix}
S_{11}&S_{22}&S_{12}&S_{13}&S_{23}
\end{bmatrix}^{T}
=\mathbf C_5\overline{\mathbf e}.}
\]

\(\mathbf C_T=\mathbf C_5\) is constant, symmetric and positive definite for the
approved material range. This is a St. Venant–Kirchhoff-type shell restriction:
material response is linear in Green–Lagrange strain and is not evidence for large
material strain validity. There are no constitutive updates or state variables.

## 10. Objective Numerical Drilling Contract

### 10.1 Element corotation

At the element center, form

\[
\mathbf J_0^c=
\begin{bmatrix}\mathbf G_\xi&\mathbf G_\eta&\mathbf G_\zeta\end{bmatrix}_{(0,0,0)},
\qquad
\mathbf j^c=
\begin{bmatrix}\mathbf g_\xi&\mathbf g_\eta&\mathbf g_\zeta\end{bmatrix}_{(0,0,0)},
\]

\[
\mathbf F_c=\mathbf j^c(\mathbf J_0^c)^{-1}.
\]

Require finite \(\det\mathbf F_c>0\). Its right polar decomposition is

\[
\mathbf F_c=\mathbf R_e\mathbf U_e,\qquad
\mathbf R_e\in SO(3),\quad
\mathbf U_e=(\mathbf F_c^T\mathbf F_c)^{1/2}\succ0 .
\]

For a superposed rigid motion,
\(\mathbf F_c^+=\mathbf R\mathbf F_c\) and
\(\mathbf R_e^+=\mathbf R\mathbf R_e\).

### 10.2 Objective drill angle

The nodal total triad rotation is

\[
\mathbf R_I=\mathbf Q_I(\mathbf Q_I^0)^T .
\]

Remove the element corotation:

\[
\mathbf C_{eI}=\mathbf R_e^T\mathbf R_I.
\]

Let \(\mathbf y_{eI}=\mathbf C_{eI}\mathbf d_I^0\). Construct the minimal swing
\(\mathbf S_{eI}\) that maps \(\mathbf d_I^0\) to \(\mathbf y_{eI}\) using the
same swing formula as §7.2, and define

\[
\mathbf T_{eI}=\mathbf S_{eI}^T\mathbf C_{eI}.
\]

\(\mathbf T_{eI}\mathbf d_I^0=\mathbf d_I^0\), so it is a twist about the initial
director. With \(\operatorname{vex}([\mathbf v]_\times)=\mathbf v\),

\[
s_{\gamma,I}=
\frac12\mathbf d_I^0\cdot
\operatorname{vex}(\mathbf T_{eI}-\mathbf T_{eI}^T),
\qquad
c_{\gamma,I}=\frac{\operatorname{tr}\mathbf T_{eI}-1}{2},
\]

\[
\boxed{
\gamma_{eI}=\operatorname{atan2}(s_{\gamma,I},c_{\gamma,I}).}
\]

The vector
\(\boldsymbol\gamma_e=[\gamma_{e1},\ldots,\gamma_{e4}]^T\in\mathbb R^4\)
contains numerical coordinates only. Antipodal swing or the non-unique twist
branch at \(\pm\pi\) fails closed.

For
\(\mathbf x=\mathbf R\mathbf X+\mathbf c\) and
\(\mathbf Q_I=\mathbf R\mathbf Q_I^0\),
\(\mathbf R_e=\mathbf R\), \(\mathbf C_{eI}=\mathbf I\), and
\(\boldsymbol\gamma_e=\mathbf0\). Therefore finite rigid motion creates no
drilling energy.

### 10.3 Fixed scale and potential

Compute the approved linear physical
\(\mathbf K_{20}^0\) once from the initial reference element. Let
\(\mathcal R_+\) contain only its eight finite positive director-tangent rotational
diagonals. Require \(\mathcal R_+\ne\varnothing\), then fix

\[
\boxed{
k_{ref,e}=\min\mathcal R_+,\qquad
k_{d,e}=10^{-3}k_{ref,e}.}
\]

\(k_{d,e}\) has force-length units and remains immutable for all increments and
iterations. It is not recomputed from a trial or committed tangent, so its residual
and Hessian have no hidden scale derivative.

\[
\boxed{
\Pi_{drill}^e=\frac12 k_{d,e}
\boldsymbol\gamma_e^T\boldsymbol\gamma_e.}
\]

It is not physical shell strain energy and creates no physical strain, stress,
resultant, or drilling output.

### 10.4 Exact drill residual and Hessian

Define

\[
\mathbf D_\gamma=
D\boldsymbol\gamma_e\in\mathbb R^{4\times24},
\qquad
\mathbf H_I^\gamma=D^2\gamma_{eI}\in\mathbb R^{24\times24}.
\]

Then

\[
\boxed{
\mathbf r_{drill}^e
=k_{d,e}\mathbf D_\gamma^T\boldsymbol\gamma_e,}
\]

\[
\boxed{
\mathbf K_{drill}^e
=k_{d,e}\left(
\mathbf D_\gamma^T\mathbf D_\gamma+
\sum_{I=1}^4\gamma_{eI}\mathbf H_I^\gamma\right).}
\]

Both derivatives include \(\mathbf R_e\), nodal exponential maps, swing removal and
atan2. A frozen-frame or secant drill tangent is prohibited.

For exact polar derivatives, write
\(\delta\mathbf F=\mathbf R_e\mathbf A\),
\(\delta\mathbf R_e=\mathbf R_e\boldsymbol\Omega\), where
\(\boldsymbol\Omega^T=-\boldsymbol\Omega\). Solve

\[
\boxed{
\boldsymbol\Omega\mathbf U_e+\mathbf U_e\boldsymbol\Omega
=\mathbf A-\mathbf A^T,}
\]

and set
\(\delta\mathbf U_e=\mathbf A-\boldsymbol\Omega\mathbf U_e\).
For mixed directions \(p,s\), define

\[
\mathbf A_p=\mathbf R_e^T\mathbf F_{c,p},\qquad
\mathbf A_{p,s}=-\boldsymbol\Omega_s\mathbf A_p+
\mathbf R_e^T\mathbf F_{c,ps},
\]

and differentiate this Sylvester equation:

\[
\boldsymbol\Omega_{ps}\mathbf U_e+
\mathbf U_e\boldsymbol\Omega_{ps}
=\mathbf A_{p,s}-\mathbf A_{p,s}^T
-\boldsymbol\Omega_p\mathbf U_{e,s}
-\mathbf U_{e,s}\boldsymbol\Omega_p ,
\]

\[
\mathbf R_{e,ps}
=\mathbf R_e(\boldsymbol\Omega_s\boldsymbol\Omega_p+
\boldsymbol\Omega_{ps}).
\]

The exchange \(p\leftrightarrow s\) must agree. Applying the norm, swing and atan2
derivative primitives from §7.3 completes exact
\(\mathbf D_\gamma,\mathbf H_I^\gamma\).

At the initial state, a pure nodal drill direction
\(\delta\boldsymbol\phi_I=\gamma_I\mathbf d_I^0\), with zero translations and
zero physical swing, gives the approved \(k_{d,e}\mathbf I_4\) action. General
deformation also contains the objective corotation coupling required to preserve
the six physical rigid modes.

## 11. Element Residual and Same-Trial Tangent

### 11.1 Physical 20-coordinate equations

\[
\Pi_{phys}^e
=\frac12\int_{\Omega_0^e}
\overline{\mathbf e}^T\mathbf C_5\overline{\mathbf e}\,dV_0.
\]

The physical residual and tangent are

\[
\boxed{
r_{20,A}^{phys}
=\int_{\Omega_0^e}\mathbf B_A^T\mathbf s\,dV_0,}
\]

\[
\boxed{
K_{20,AB}^{mat}
=\int_{\Omega_0^e}
\mathbf B_A^T\mathbf C_5\mathbf B_B\,dV_0,}
\]

\[
\boxed{
K_{20,AB}^{geo}
=\int_{\Omega_0^e}
\mathbf s^T\mathbf G_{AB}\,dV_0,}
\]

\[
\mathbf r_{20}^{phys}\in\mathbb R^{20},\qquad
\mathbf K_{20}^{phys}
=\mathbf K_{20}^{mat}+\mathbf K_{20}^{geo}
\in\mathbb R^{20\times20}.
\]

The geometric term includes the exact director Hessian and the quadratic
Green–Lagrange terms. It is not a separately guessed initial-stress matrix.

### 11.2 Fixed-spatial nodal load in physical coordinates

In the element equations below,
\(\mathbf F_I^\star,\mathbf M_{I,src}^\star\) denote a deterministic
source-ordered element load share whose incident-element sum is the single
aggregated nodal CLOAD. Equivalently the identical nodal contribution may be
assembled once outside element kernels. It must never be copied in full to every
incident element. Define the per-unit-load-factor effective physical covector

\[
\mathbf M_{I,phys}^\star(\mathbf q)
:=\mathbf P_I(\mathbf q)\mathbf M_{I,src}^\star,
\qquad
\mathbf M_{I,phys}(\mathbf q;\lambda)
=\lambda\mathbf M_{I,phys}^\star(\mathbf q).
\]

For nodal physical swing coordinates
\(\chi_{I,a}^2,\ a\in\{1,2\}\), define

\[
\mathbf L_{I,a}
=\mathbf d_I\times
\frac{\partial\mathbf d_I}{\partial\chi_{I,a}^2}
\in\mathbb R^3,
\qquad
m_{I,a}^{20}
=\mathbf M_{I,phys}^\star\cdot\mathbf L_{I,a}
=\mathbf M_{I,src}^\star\cdot\mathbf L_{I,a}.
\]

Insert these two entries per node, and zeros in translation positions, into
\(\mathbf m_{20}^{M}(\mathbf q_{20})\in\mathbb R^{20}\).
Since \(\mathbf d_I\cdot\mathbf L_{I,a}=0\), the two displayed forms of
\(m_{I,a}^{20}\) are exactly identical even when the admitted source has a small
director component. Its exact physical-chart derivative is

\[
\begin{aligned}
K_{I,ab}^{M,20}
&=\frac{\partial\mathbf M_{I,phys}^\star}
{\partial\chi_{I,b}^2}\cdot\mathbf L_{I,a}
+\mathbf M_{I,phys}^\star\cdot
\frac{\partial\mathbf L_{I,a}}{\partial\chi_{I,b}^2}\\
&=\mathbf M_{I,src}^\star\cdot\left(
\frac{\partial\mathbf d_I}{\partial\chi_{I,b}^2}
\times
\frac{\partial\mathbf d_I}{\partial\chi_{I,a}^2}
+
\mathbf d_I\times
\frac{\partial^2\mathbf d_I}
{\partial\chi_{I,a}^2\partial\chi_{I,b}^2}
\right).
\end{aligned}
\]

The first line explicitly contains
\(\partial\mathbf M_{I,phys}^\star/\partial\chi_{I,b}^2
=(\partial\mathbf P_I/\partial\chi_{I,b}^2)
\mathbf M_{I,src}^\star\). The second line is the algebraically identical
immutable-source form obtained by differentiating
\(\mathbf M_{I,src}^\star\cdot\mathbf L_{I,a}\); therefore the projector derivative
must not be dropped when the first form is used.

All cross-node and translation derivatives are zero. In general
\(\mathbf K_{20}^{M}\) need not be symmetric; this is the exact derivative of a
constant-spatial source-moment work one-form restricted to physical swing space.

Let
\(\mathbf f_{24}^{F,\star}\in\mathbb R^{24}\) be the translational dead-force
selection vector. The element external generalized vector and its exact derivative
are

\[
\boxed{
\mathbf f_{ext}^e
=\lambda\left(
\mathbf f_{24}^{F,\star}
+\mathbf A^T\mathbf m_{20}^{M}\right),}
\]

\[
\boxed{
\mathbf K_{ext}^{M,e}
=\lambda\left[
\mathbf A^T\mathbf K_{20}^{M}\mathbf A
+\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi
\right].}
\]

As an independent exact global-coordinate identity, the per-unit-load-factor nodal
rotational block is

\[
\mathbf f_{\phi,I}^{M,\star}
=\mathbf J_l(\boldsymbol\phi_I)^T
\mathbf P_I\mathbf M_{I,src}^\star
=\mathbf J_l(\boldsymbol\phi_I)^T\mathbf M_{I,phys}^\star,
\qquad
\mathbf P_I=\mathbf I-\mathbf d_I\mathbf d_I^T,
\]

\[
D\mathbf f_{\phi,I}^{M,\star}[\mathbf p]
=D\mathbf J_l(\boldsymbol\phi_I)[\mathbf p]^T
\mathbf P_I\mathbf M_{I,src}^\star
+\mathbf J_l(\boldsymbol\phi_I)^T
D\mathbf P_I[\mathbf p]\mathbf M_{I,src}^\star,
\]

\[
D\mathbf P_I[\mathbf p]
=-\left(
D\mathbf d_I[\mathbf p]\mathbf d_I^T+
\mathbf d_I D\mathbf d_I[\mathbf p]^T\right).
\]

Let \(\mathbf f_{24}^{M,\star}\) contain these four nodal rotational blocks and
zero translation blocks. Virtual-work covariance of the exact map gives

\[
\boxed{
\mathbf A^T\mathbf m_{20}^{M}=\mathbf f_{24}^{M,\star},
\qquad
D(\mathbf A^T\mathbf m_{20}^{M})[\mathbf p]
=\left(\mathbf A^T\mathbf K_{20}^{M}\mathbf A
+\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi\right)\mathbf p
=D\mathbf f_{24}^{M,\star}[\mathbf p].}
\]

Thus the physical-chart and global-coordinate forms are the same residual
covector and the same derivative of the §5 work one-form. Both include the
projector derivative; neither changes the immutable source identity nor creates
follower-load direction updates or a drilling load channel.

### 11.3 Complete 24-coordinate equations

The element residual is

\[
\boxed{
\mathbf r_e=
\mathbf A^T\mathbf r_{20}^{phys}
+\mathbf r_{drill}^e
-\lambda\left(
\mathbf f_{24}^{F,\star}
+\mathbf A^T\mathbf m_{20}^{M}\right)
\in\mathbb R^{24}.}
\]

The exact same-trial algorithmic tangent is

\[
\boxed{
\begin{aligned}
\mathbf K_{T,e}
=&\ \mathbf A^T\mathbf K_{20}^{phys}\mathbf A\\
&+\sum_{a=1}^{20}r_{20,a}^{phys}\mathbf H_a^\Phi
+\mathbf K_{drill}^e\\
&-\lambda\left[
\mathbf A^T\mathbf K_{20}^{M}\mathbf A
+\sum_{a=1}^{20}m_{20,a}^{M}\mathbf H_a^\Phi
\right]
\in\mathbb R^{24\times24}.
\end{aligned}}
\]

The first line contains mapped material and Green–Lagrange geometric terms. The
second line contains physical coordinate-map curvature and the objective drill
gradient/Hessian. The last two lines are the complete fixed-spatial source-moment
work-conjugacy derivative, including the effective-work projector derivative,
rotation-coordinate derivative and map-Hessian term. Fixed-global nodal force has
no configuration derivative.

No term may be evaluated from a committed or previous-iteration state while another
term uses the current trial. In particular, the following same-trial identities
are mandatory:

\[
\mathbf K_{T,e}=D\mathbf r_e,\qquad
\mathbf K_T=D\mathbf r,\qquad
\mathbf r=\sum_e\mathcal A_e^T\mathbf r_e .
\]

### 11.4 Global equations and reactions

Assemble element contributions and nodal load shares in stable source order:

\[
\mathbf r(\mathbf q;\lambda)
=\mathbf f_{int}(\mathbf q)
-\mathbf f_{ext}(\mathbf q;\lambda),
\qquad
\mathbf K_T=D_{\mathbf q}\mathbf r.
\]

At a trial that satisfies its prescribed coordinates, Full Newton on free
coordinates is

\[
\boxed{
\mathbf K_{ff}^{(k)}\Delta\mathbf q_f^{(k)}
=-\mathbf r_f^{(k)},\qquad
\Delta\mathbf q_c^{(k)}=\mathbf0.}
\]

The tangent is freshly factorized and solved once for every entered iteration.
A valid \(0\times0\) free system is a trivial solve, not singular.

The full residual is preserved. Its constrained entries are generalized reactions
conjugate to global \((\mathbf u,\boldsymbol\phi)\); free entries are equilibrium
evidence. Translation entries are spatial forces. For balance or global physical
moment reporting, convert a rotational generalized covector by

\[
\boxed{
\mathbf m_{sp,I}
=\mathbf J_l(\boldsymbol\phi_I)^{-T}\mathbf r_{\phi,I}.}
\]

This conversion is applied consistently to internal, external and reaction
covectors; it does not alter the Newton equations.

## 12. Mapping and Numerical Integration

### 12.1 Reference Jacobian and derivative transformation

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

\[
dV_0=J_0\,d\xi\,d\eta\,d\zeta,\qquad
\nabla_0(\cdot)=(\cdot)_{,\alpha}\otimes\mathbf G^\alpha .
\]

Every derivative transformation uses the inverse/contravariant basis of this fixed
reference Jacobian. The current mapping

\[
\mathbf j=
\begin{bmatrix}\mathbf g_\xi&\mathbf g_\eta&\mathbf g_\zeta\end{bmatrix},
\qquad j=\det\mathbf j
\]

is not the integration measure but must remain finite and positive to exclude local
inversion.

### 12.2 Required validation locations

Validate finite data, nonzero surface area, finite \(J_0>0\), and finite
\(j>0\) at:

- all eight volume integration points;
- all four midsurface MITC tying points;
- the element center used by drilling corotation;
- every bottom/middle/top recovery point before commit.

Also require finite reciprocal bases, nonzero interpolated director, valid fixed
local frame, \(\det\mathbf F_c>0\), SPD polar stretch, valid swing/twist branches,
and \(Q^TQ=I,\det Q=+1\). A failed point is not omitted, averaged, clamped, or
replaced.

### 12.3 Quadrature

Use tensor-product Gauss points

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

Thus each element uses four midsurface by two thickness points, eight volume
evaluations. The following are integrated at all eight points with the same
reference \(J_0w\):

- physical internal residual;
- material tangent;
- geometric tangent;
- physical strain-energy and physical recovery thickness moments.

The MITC tying values are evaluated at the four \(\zeta=0\) tying locations and
used for both thickness points. The objective drill potential is an element-level
four-coordinate numerical term and is not multiplied by volume quadrature.

S4 and S4R use this identical rule. One-point integration, reduced/selective
integration, source-selected hourglass behavior, and higher-order production
quadrature are prohibited.

## 13. Output Recovery

Only accepted committed states are recoverable authoritative frames. No trial
Newton state is an authoritative result.

### 13.1 Nodal kinematics and reactions

At each source node:

- displacement \([U1,U2,U3]=\mathbf u_I\), global, length;
- rotation \([UR1,UR2,UR3]=\boldsymbol\phi_I\), fixed-initial global rotation-vector
  components, dimensionless;
- committed \(\mathbf Q_I\in SO(3)\) and \(\mathbf d_I=\mathbf Q_I\mathbf e_3\) as
  the kinematic state identity;
- translation reaction from constrained full residual entries, global force;
- rotational generalized reaction and its spatial conversion from §11.4,
  force-length.

For unambiguous recovery, let
\(\mathbf C_{u,I},\mathbf C_{\phi,I}\in\mathbb R^{3\times3}\) be diagonal
selectors with one on a prescribed global component and zero otherwise. With the
accepted full-residual blocks
\(\mathbf r_{u,I},\mathbf r_{\phi,I}\in\mathbb R^3\), define

\[
\boxed{
\mathbf F_I^{react}:=\mathbf C_{u,I}\mathbf r_{u,I}
\in\mathbb R^3,
\qquad
\mathbf m_{\phi,I}^{react}:=\mathbf C_{\phi,I}\mathbf r_{\phi,I}
\in\mathbb R^3,}
\]

\[
\boxed{
\mathbf M_I^{react}
:=\mathbf J_l(\boldsymbol\phi_I)^{-T}\mathbf m_{\phi,I}^{react}
\in\mathbb R^3.}
\]

\(\mathbf F_I^{react}\) has force units and is work-conjugate to
\(\delta\mathbf u_I\). The generalized covector
\(\mathbf m_{\phi,I}^{react}\) and spatial reaction moment
\(\mathbf M_I^{react}\) have force-length units and satisfy
\(\mathbf m_{\phi,I}^{react}\cdot\delta\boldsymbol\phi_I
=\mathbf M_I^{react}\cdot\delta\boldsymbol\vartheta_I\).
Unselected free residual components remain equilibrium evidence and are not
relabelled as reactions. The symbol \(\mathbf R_I\) remains reserved exclusively
for the rotation matrix defined in §§3.2 and 10.2.

Exact HDF5 projection is I/O Definition ownership. This section fixes the quantity,
coordinate and unit meaning only.

### 13.2 Integration-point physical strain and stress

At each four midsurface Gauss locations and both thickness Gauss points recover
\(\overline{\mathbf e}\) and \(\mathbf s\) from the exact accepted state, same
MITC tying, fixed reference frame and material used by residual/tangent.

Generalized strains are thickness moments:

\[
\boldsymbol\epsilon_0
=\frac1t\int_{-t/2}^{t/2}
\begin{bmatrix}E_{11}&E_{22}&2E_{12}\end{bmatrix}^T dz,
\]

\[
\boldsymbol\kappa
=\frac{12}{t^3}\int_{-t/2}^{t/2}
z\begin{bmatrix}E_{11}&E_{22}&2E_{12}\end{bmatrix}^T dz,
\]

\[
\boldsymbol\gamma_0
=\frac1t\int_{-t/2}^{t/2}
\begin{bmatrix}2E_{13}&2E_{23}\end{bmatrix}^T dz .
\]

Their order is

\[
[E11,E22,G12,K11,K22,K12,G13,G23].
\]

Membrane/shear strain is dimensionless and curvature has \(1/\text{length}\).

### 13.3 Section resultants and stresses

\[
\mathbf N=\int_{-t/2}^{t/2}
\begin{bmatrix}S_{11}&S_{22}&S_{12}\end{bmatrix}^T dz,
\]

\[
\mathbf M=\int_{-t/2}^{t/2}
z\begin{bmatrix}S_{11}&S_{22}&S_{12}\end{bmatrix}^T dz,
\qquad
\mathbf Q_s=\int_{-t/2}^{t/2}
\begin{bmatrix}S_{13}&S_{23}\end{bmatrix}^T dz .
\]

The order remains

\[
[N11,N22,N12,M11,M22,M12,Q13,Q23].
\]

\(\mathbf N,\mathbf Q_s\) have force/length units and \(\mathbf M\) has force
units. Positive \(z\) and positive bending sign are the approved linear convention.
At \(z=-t/2,0,+t/2\), recover only \([S11,S22,S12]\) as direct evaluations.
\(S_{33}=0\) is documented but not emitted.

There is no nodal extrapolation/averaging and no numerical drill stress, resultant,
energy ratio, or output. Exact element, natural-coordinate, thickness-position,
frame and source identity remain attached to every physical row.

### 13.4 Equilibrium identities

At an accepted state,

\[
\mathbf r_f=\mathbf0
\]

in the exact equation; numerical acceptance tolerance is downstream. Define the
applied physical load resultants at the accepted state by

\[
\mathbf F_I^{appl}(\lambda):=\lambda\mathbf F_I^\star
\in\mathbb R^3,
\qquad
\mathbf M_I^{appl}(\mathbf q;\lambda)
:=\mathbf M_{I,phys}(\mathbf q;\lambda)
=\lambda\mathbf P_I\mathbf M_{I,src}^\star
\in\mathbb R^3.
\]

Here \(\mathbf F_I^{appl}\) has force units and
\(\mathbf M_I^{appl}\) has force-length units. For a fixed global origin
\(\mathbf X_o\), physical resultant evidence uses current nodal positions
\(\mathbf x_I=\mathbf X_I+\mathbf u_I\):

\[
\boxed{
\mathbf e_F
=\sum_I(\mathbf F_I^{appl}+\mathbf F_I^{react})
\in\mathbb R^3\ \text{[force]},}
\]

\[
\boxed{
\mathbf e_M=\sum_I[
(\mathbf x_I-\mathbf X_o)\times
(\mathbf F_I^{appl}+\mathbf F_I^{react})
+\mathbf M_I^{appl}+\mathbf M_I^{react}]
\in\mathbb R^3\ \text{[force-length]}.}
\]

The balance uses the exact effective physical-work moment, while the immutable
source \(\mathbf M_{I,src}^\star\) remains unchanged for source identity and the
fixed-global load path. Its removed director component is neither a physical
balance contribution nor a numerical drill load.
The objective drilling potential contributes self-equilibrated nodal forces/moments
and is included in full residual reactions; it has no separate physical output.
Balance acceptance tolerance is owned by Numerical Review.

## 14. Newton and State Algorithm Pseudocode

~~~text
preprocess in stable source order:
    validate source S4/S4R mapping, material, section, node order and initial geometry
    construct approved initial nodal triads Q0 and directors d0
    aggregate immutable fixed-global nodal CLOAD sources and validate initial moment rho_M
    validate all reference Gauss/tying/center/recovery mappings
    compute initial linear physical K20_0 for each element
    select k_ref from finite positive physical tangent-rotation diagonals
    freeze k_d = 1e-3 * k_ref for the entire analysis

initialize:
    committed = exact step-start state at factor 0
    intended factors follow the approved fixed increment and final shortening

for each intended factor lambda in deterministic order:
    trial = deep copy(committed)
    set every constrained translation and rotation coordinate to its absolute path value
    keep the preceding committed free coordinates as the initial free guess

    for iteration k = 1..100:
        derive every Q = Exp(phi^) * Q0 and d = Q * e3 from this trial
        validate finite state, principal charts, current geometry and source-moment rho_M

        for each element in stable source order:
            evaluate one trial identity at all Gauss/tying/center locations
            form MITC Green-Lagrange strain and exact first/second variations
            form r20_phys, K20_mat and K20_geo
            form exact Phi Jacobian A and all map Hessians H_Phi
            form polar corotation, objective gamma, D_gamma and H_gamma
            form r_drill and exact K_drill
            form fixed-force and projected physical-moment work plus exact P/J_l/A derivatives
            form complete r_e and K_T_e from §11.3

        deterministically assemble full r and K_T
        restrict to the stable free coordinate order
        freshly factorize K_ff for this iteration
        solve K_ff * Delta_q_f = -r_f; set Delta_q_c = 0
        reject any nonfinite residual, tangent, factorization or correction

        evaluate all four inclusive convergence gates from:
            current same-trial free residual, internal force and external force
            fresh same-iteration free correction
            current trial incremental solution relative to step-start committed state
        emit exactly one deterministic diagnostic summary for entered iteration k

        if all four gates pass:
            recover all physical rows from this unchanged current trial
            validate complete finite candidate and state/output identities
            atomically commit displacement, phi, Q, d, factor, forces, reaction and rows
            create exactly one accepted-frame candidate identity
            exit the iteration loop

        if k == 100:
            discard trial and fail without iteration 101, cutback or retry

        apply the complete fresh correction unchanged:
            u <- u + Delta_u
            phi <- phi + Delta_phi
            Q_new <- Exp((phi_new)^) * Q0 = Delta_R * Q_old
        continue with a new same-trial evaluation

on any increment/recovery/output failure:
    discard the complete trial/candidate
    preserve prior committed state and prior authoritative final result
    perform no replacement increment or hidden solver fallback
~~~

### 14.1 Four convergence channels

Let free translation and rotation-coordinate subvectors carry superscripts \(u\)
and \(\phi\). Residual channels use

\[
\eta_R^u=
\frac{\|\mathbf r_f^u\|_\infty}
{\max(\|\mathbf f_{int,f}^u\|_\infty,
\|\mathbf f_{ext,f}^u\|_\infty)},
\]

\[
\eta_R^\phi=
\frac{\|\mathbf r_f^\phi\|_\infty}
{\max(\|\mathbf f_{int,f}^\phi\|_\infty,
\|\mathbf f_{ext,f}^\phi\|_\infty)}
\]

when the matching denominator is nonzero. They pass iff
\(\eta_R^u\le10^{-2}\) and \(\eta_R^\phi\le10^{-2}\).
An exact-zero denominator instead passes iff the numerator is at most
\(10^{-3}\) model-consistent force or \(10^{-3}\) model-consistent force-length.

Let

\[
\Delta\mathbf q_{inc}
=\mathbf q_{trial}-\mathbf q_{step-start}.
\]

Correction channels use

\[
\eta_\Delta^u=
\frac{\|\Delta\mathbf u_f\|_\infty}
{\|\Delta\mathbf u_{inc,f}\|_\infty},
\qquad
\eta_\Delta^\phi=
\frac{\|\Delta\boldsymbol\phi_f\|_\infty}
{\|\Delta\boldsymbol\phi_{inc,f}\|_\infty}
\]

for nonzero denominators and pass iff each is \(\le10^{-2}\). An exact-zero
denominator instead passes iff the fresh numerator is at most \(10^{-3}\) length
or \(10^{-3}\) dimensionless rotation. Empty free channels have exact-zero
numerator/denominator and pass their typed zero-scale rule.

All four gates are combined by logical AND. Constrained coordinates are excluded.
These are approved solver-convergence values, not external-reference tolerances.

## 15. Verification Identities

The following are formulation identities whose numerical norms/tolerances are set
by Numerical Review.

### 15.1 Linear limit

At \(\mathbf q=\mathbf0\), stress-free reference:

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

\[
\mathbf K_{phys}^{24}
=\mathbf T_p^T\mathbf K_{20}^{linear}\mathbf T_p.
\]

The nonlinear Green–Lagrange/MITC physical strain, stress, resultants and recovery
reduce to the approved linear MITC4 signs, component orders and 2 x 2 x 2 rule.
A pure drill basis reduces to \(k_d\mathbf I_4\). Objective corotation adds only
the coupling required to make finite rigid motion drill-free and contributes
neither physical recovery nor stress.

### 15.2 Objectivity and six rigid modes

For any \(\mathbf R\in SO(3)\), \(\mathbf c\in\mathbb R^3\),

\[
\mathbf x^+=\mathbf R\mathbf x+\mathbf c,\qquad
\mathbf Q_I^+=\mathbf R\mathbf Q_I
\]

implies

\[
\overline{\mathbf E}^+=\overline{\mathbf E},\quad
\Pi_{phys}^+=\Pi_{phys},\quad
\boldsymbol\gamma_e^+=\boldsymbol\gamma_e,\quad
\Pi_{drill}^+=\Pi_{drill}.
\]

At a stress-free reference with no loads, three translations and three proper
rigid rotations produce zero residual and zero energy. After objective drill
regularization, an isolated valid element has exactly these six physical rigid
modes; pure numerical drill is not an additional zero mode.

### 15.3 Tangent consistency

For a finite valid trial and deterministic direction \(\mathbf p\),

\[
\boxed{
\mathbf K_T(\mathbf q)\mathbf p
=D\mathbf r(\mathbf q)[\mathbf p].}
\]

The directional comparison must exercise separately and together:

- material term;
- Green–Lagrange geometric term;
- \(\mathbf A^T\mathbf K_{20}\mathbf A\);
- \(\sum r_{20,a}\mathbf H_a^\Phi\);
- fixed-spatial source-moment projector, rotation-coordinate and map-Hessian terms;
- objective drill gradient and exact Hessian.

The two perturbed states share the same fixed reference, load factor and committed
base. Perturbation schedule and tolerance are Numerical Review ownership.

### 15.4 Symmetry, equilibrium and state

- With no applied nodal moment, the complete unreduced physical/drill tangent is
  symmetric in exact arithmetic. With a fixed-spatial nodal moment, the exact
  dead-moment work tangent is generally nonsymmetric and must not be symmetrized.
- Element and assembled internal actions satisfy virtual-work covariance under
  stable scatter/assembly.
- Accepted free residual is equilibrium evidence; constrained full residual is
  reaction.
- Global resultant force and moment identities in §13.4 include nodal couples and
  current force arms.
- Every committed factor has one accepted frame; no failed trial changes the
  preceding \(\mathbf u,\boldsymbol\phi,\mathbf Q,\mathbf d\), force, reaction or
  recovery identity.

### 15.5 Patch and recovery checks

- partition of unity, Kronecker delta and derivative-sum identities;
- constant membrane, bending, shear and twist continuation from the linear limit;
- direct/tied shear component identity and sign;
- bottom/middle/top stress and section-resultant thickness-moment checks;
- regular, rigidly rotated and admissibly distorted patches;
- no physical strain/stress/resultant change from pure numerical drill;
- stable source, element, Gauss, section and component ordering.

## 16. Numerical Risks

| risk | consequence | formulation control / downstream evidence |
| --- | --- | --- |
| Rigid-body mode loss | artificial force/energy | objective physical and drill potentials; exactly six physical rigid modes |
| Patch-test failure | inconsistent interpolation/sign | direct membrane/bending plus identical MITC projection of strain and two variations |
| Tangent asymmetry/inconsistency | degraded or failed Full Newton | exact material, geometric, map, drill Hessians and exact dead-moment work derivative; directional derivative review |
| Indefinite tangent | valid near-buckling state or numerical failure | do not require positive definiteness globally; distinguish expected physical indefiniteness from singular/malformed tangent |
| Hourglass / zero-energy mode | singular element | no reduced integration; 2 x 2 x 2 full rule; objective drill removes only numerical drill modes |
| Transverse-shear locking | overly stiff thin response | canonical MITC tying; thickness/mesh characterization |
| Membrane locking | overly stiff curved/distorted original MITC4 | documented known limitation; do not claim MITC4+ performance |
| Volumetric locking | not applicable to approved plane-stress shell | no 3-D incompressible claim; material range remains \(\nu<0.5\) |
| Director/swing antipode | non-unique physical chart | fail closed at \(d^0\cdot d=-1\); no silent axis choice |
| Exponential/log branch | non-unique global UR target/state | principal \(\|\phi\|<\pi\); downstream branch tolerance |
| Drill twist branch | discontinuous numerical potential | fail closed at non-unique \(\pm\pi\) twist |
| Nonobjective drilling | artificial rigid energy | polar element corotation and exact potential Hessian |
| State-dependent drill scale | missing scale derivative/nondeterminism | freeze initial approved \(k_d\); no trial recomputation |
| Director-parallel applied moment | physical work mixed with drill | immutable-source current-trial fail-closed ratio; exact physical-work projection; no drill channel |
| Moment derivative omission | inconsistent residual/tangent | immutable \(M_{src}^\star\), exact projector/coordinate/map derivative, no follower direction |
| Reversed/singular Jacobian | invalid basis, sign or inverted shell | \(J_0>0\), \(j>0\), \(\det F_c>0\) at every required location |
| Distortion/warpage | conditioning and membrane locking | exact point inventory; no omitted/averaged failed location |
| Nonfinite trial | corrupted state/output | reject before solve acceptance or commit |
| Constraint/target confusion | wrong load path or reaction | absolute path applied before iteration; constrained Newton correction exactly zero |
| Mixed units in convergence | false pass/fail | force/moment and translation/rotation channels remain separate |
| Determinism loss | irreproducible reference/evidence | stable source order, local buffers, fixed reduction and diagnostic order |

## 17. Requirement and Research Traceability

| requirement / gap | formulation closure |
| --- | --- |
| FESA-REQ-NSNR-003 | §§2, 6–10: original-style MITC4, homogeneous isotropic materially linear TL shell |
| FESA-REQ-NSNR-041 | §§2, 6–9, 12: one fixed initial reference, Green–Lagrange / second-Piola pair |
| FESA-REQ-NSNR-074 | §3: 6 global DOF, dimensionless correction, committed/trial \(Q,d\), multiplicative update |
| FESA-REQ-NSNR-006 | §§11, 14: complete same-trial residual/tangent and fresh Full Newton solve |
| FESA-REQ-NSNR-007/069/070 | §14.1: four typed inclusive gates and constrained exclusion |
| FESA-REQ-NSNR-008/014 | §14: finite data, iterations 1..100, unchanged correction, no fallback |
| FESA-REQ-NSNR-046/048 | §§4, 14: one proportional factor and fixed-step lifecycle |
| FESA-REQ-NSNR-018 | §§4.4, 5, 11.2: immutable fixed-global source coefficients, exact tangent-space work covector/derivative and no follower direction |
| FESA-REQ-NSNR-063 | §§4.2–4.3: factor-one absolute target and zero constrained Newton correction |
| FESA-REQ-NSNR-011 | §§3.5, 14: separate trial and atomic complete commit/rollback |
| FESA-REQ-NSNR-019 | §§11.4, 13.1, 13.4, 14: deterministic full residual, dimensioned reaction covectors and force/moment balance |
| FESA-REQ-NSNR-030 | §15: linear limit, tangent, equilibrium, objectivity and rollback identities |
| O-ROT-01 | §§3.2–3.4: fixed-reference exponential coordinate; exact left-spatial multiplicative update and derivatives |
| O-ROT-02 | §§3, 6–8, 12: nodal triad/director, exact director interpolation, point frame/branch validation |
| O-MAP-01 | §§7, 11: exact \(\Phi:\mathbb R^{24}\to\mathbb R^{20}\), \(A\), \(H^\Phi\), complete map curvature |
| O-DRILL-01 | §10: objective polar-corotated twist potential, frozen scale, exact residual/Hessian, six rigid modes |
| O-DRILL-02 | §§4.4, 5, 10, 11.2: current-director source admissibility, exact physical-work projection and no drill load channel |
| O-MOM-01 | §§5, 11.2–11.3: immutable fixed-spatial source, projected physical swing covector, exact work derivative and no follower direction |
| O-BC-01 | §§3.2–3.3, 4.2–4.3, 14: absolute component path, exact multiplicative triad update, \(\Delta q_c=0\) |

Theory basis follows research sources S1/S6 for fixed-reference nonlinear
continuum-based MITC shell and tying/integration, S2–S5 for finite rotation manifold
and exact potential derivatives, and S7 for same-state Full Newton. Local wiki
cross-checks were [[Total Lagrangian Shell Formulation]],
[[Green-Lagrange Strain Linearization]], [[MITC4 Shell Element]],
[[MITC Shell Kinematics]], [[Continuum Mechanics Based Four-Node Shell Element]],
[[Assumed Transverse Shear Strain Interpolation]], and
[[Nonlinear Finite Element Analysis]]. Project requirements and approved repository
contracts take precedence over those navigation pages.

## 18. Open Issues and Downstream Handoff

### 18.1 Formulation closures

O-ROT-01, O-ROT-02, O-MAP-01, O-DRILL-01, O-DRILL-02, O-MOM-01 and O-BC-01
are closed by the explicit equations and fail-closed boundaries in this document.
No placeholder derivation remains for implementation of their residual or tangent.

### 18.2 Numerical Review-owned decisions

The following remain intentionally open to Numerical Review and do not alter this
formulation:

- O-TOL-01: scale-aware tolerances for tangent derivative, symmetry, objectivity,
  triad invariants, equilibrium, branch proximity and Jacobian validity;
- O-REF-01/O-REF-02: exact stored nonlinear input/CSV and independent-source
  readiness;
- O-BENCH-01: exact blocking large-response benchmark and quantities;
- O-STEP-01: endpoint arithmetic and fixed increment adequacy without cutback;
- O-MMS-01: whether shell MMS remains nonblocking.

Numerical Review must independently check dimensions, signs, exact derivative
identities, the polar-corotation drill construction, six rigid modes, prescribed
rotation path, moment rejection, tangent symmetry/consistency, locking limitations,
and benchmark traversability. This status is ready-for-numerical-review, not
numerically approved and not a downstream implementation authorization.

### 18.3 Boundary of later stages

- Exact HDF5 datasets, log representation, parser/semantic validation and public CLI
  behavior remain I/O Definition ownership.
- Exact reference artifacts, comparison families and tolerance application remain
  Numerical Review/reference contract ownership under docs/TOLERANCE.md.
- C++ type/API/file ownership, sparse implementation and tests remain later planning
  and implementation ownership.
- No Abaqus, Nastran or other reference solver execution, artifact creation,
  artifact modification, or internal-equivalence claim is authorized by this
  document.