FESADev/docs/MITC4_FORMULATION.md

# MITC4 Formulation

## Purpose
This document is the formulation contract for FESA's Phase 1 four-node MITC4 shell element.

It supersedes the earlier simplified Phase 1 baseline that used an averaged-edge shell basis, analytic thickness integration, and `drilling_stiffness_scale = 1.0e-6`. The next Phase 1 implementation pass must treat this document as the MITC4 gate before coding or evaluating shell stiffness.

The goal is not to reproduce Abaqus internals. The goal is to implement a documented, testable Dvorkin-Bathe-style MITC4 element that can be compared against stored Abaqus `S4` reference artifacts without running Abaqus locally.

## Local Source Basis
The formulation below is based on the local papers in `docs/Paper/`:

- `docs/Paper/AContinuumMechanicsBasedFourNodeShellElementforGeneralNonlinearAnalysis/`
  - Dvorkin and Bathe's original four-node non-flat shell element.
  - Defines the degenerated-continuum geometry, displacement interpolation, convected transverse shear interpolation, constitutive transformation, and benchmark expectations.
- `docs/Paper/FourNodeQuadrilateralShellElementMITC4/`
  - Restates the MITC4 geometry and displacement interpolation, gives explicit midside shear tying relations, and reports patch-test and Scordelis-Lo verification behavior.
- `docs/Paper/유한요소해석법을이용한쉘구조물의동적좌굴해석/`
  - Provides the Korean thesis derivation for covariant strain, plane-stress material matrix with shear correction, 6-DOF shell transformation, drilling stabilization, and nonlinear extension notes.
- `docs/Paper/2007쉘구조물의유한요소해석에대하여/`
  - Provides shell-model background, locking behavior, asymptotic behavior, and benchmark/test interpretation.
- `docs/Paper/mitc공부/`
  - Treated as supporting study notes only. Use the papers above for binding formulas.

## Phase 1 Scope
Phase 1 implements:

- Abaqus `TYPE=S4` mapped to FESA `MITC4`.
- Four-node quadrilateral shell in Abaqus S4 node order.
- Linear static, small-strain, linear isotropic elastic analysis.
- Homogeneous shell section with one thickness value per element for the initial implementation.
- Six global DOFs per node: `UX`, `UY`, `UZ`, `RX`, `RY`, `RZ`.
- MITC transverse shear interpolation in a convected coordinate system.
- Constrained DOF elimination, full-vector reaction recovery, and minimum output fields `U` and `RF`.

Phase 1 does not implement:

- Abaqus `S4R`, reduced integration, or hourglass control.
- Geometric nonlinearity, material nonlinearity, dynamics, pressure loads, thermal-stress coupling, composite sections, or mesh quality diagnostics.
- Abaqus `*Orientation`, user-defined nodal normals, or shell offset behavior unless added by a later ADR.

## Natural Coordinates and Node Order
FESA uses the standard bilinear quadrilateral natural coordinates:

```text
node 1: (xi, eta) = (-1, -1)
node 2: (xi, eta) = (+1, -1)
node 3: (xi, eta) = (+1, +1)
node 4: (xi, eta) = (-1, +1)
zeta in [-1, +1] through the thickness
```

Shape functions:

```text
N1 = 0.25 * (1 - xi) * (1 - eta)
N2 = 0.25 * (1 + xi) * (1 - eta)
N3 = 0.25 * (1 + xi) * (1 + eta)
N4 = 0.25 * (1 - xi) * (1 + eta)
```

Midside MITC tying points in the FESA convention:

```text
A = ( 0, -1)  edge 1-2
B = (-1,  0)  edge 1-4
C = ( 0, +1)  edge 4-3
D = (+1,  0)  edge 2-3
```

Dvorkin-Bathe papers may use a natural-coordinate orientation where the signs on the `A/C` interpolation are reversed relative to the FESA convention above. FESA implementation must preserve the edge definitions above and use the FESA interpolation formulas in the MITC shear section.

## Degenerated-Continuum Geometry
The shell is represented as a 3D continuum degenerated through the thickness:

```text
X(xi, eta, zeta) =
  sum_k N_k(xi, eta) X_k
  + zeta / 2 * sum_k t_k N_k(xi, eta) Vn_k
```

where:

- `X_k` is the midsurface coordinate of node `k`.
- `t_k` is the nodal thickness. Phase 1 uses `t_k = t` for all four nodes.
- `Vn_k` is the nodal director. The papers allow the director to be non-normal to the midsurface; Phase 1 derives it from geometry because Abaqus nodal normals are not in the parser subset.

Phase 1 director policy:

- If no parsed nodal directors exist, use the element-center midsurface normal for all four nodal directors:

```text
G1_c = dX_mid / dxi  at (0, 0)
G2_c = dX_mid / deta at (0, 0)
Vn_k = normalize(G1_c x G2_c), k = 1..4
```

- If `G1_c x G2_c` is near zero, reject the element with an invalid/singular element diagnostic.
- This flat-or-mildly-warped policy is sufficient for the Phase 1 S4 baseline. A true non-flat nodal-director policy requires a later ADR and benchmark update.

The nodal local director axes are:

```text
V1_k = normalize(EY x Vn_k)
V2_k = Vn_k x V1_k
```

If `EY x Vn_k` is near zero, use a deterministic fallback axis before normalization:

```text
V1_k = normalize(EZ x Vn_k)
if still near zero, V1_k = normalize(EX x Vn_k)
V2_k = Vn_k x V1_k
```

The implementation must keep `V1_k`, `V2_k`, and `Vn_k` orthonormal and right-handed.

## Displacement and Rotation Interpolation
The local five-DOF MITC4 displacement field is:

```text
u(xi, eta, zeta) =
  sum_k N_k u_k
  + zeta / 2 * sum_k t_k N_k q_k

q_k = -V2_k * alpha_k + V1_k * beta_k
```

where:

- `u_k` is the nodal translation vector.
- `alpha_k` is the rotation about `V1_k`.
- `beta_k` is the rotation about `V2_k`.
- Drilling rotation `gamma_k` about `Vn_k` does not enter the continuum strain field.

FESA stores six global DOFs per node. Global nodal rotations are mapped to local director rotations by:

```text
theta_k = [RX_k, RY_k, RZ_k]^T

[alpha_k, beta_k, gamma_k]^T = L_k theta_k

L_k =
  [ EX dot V1_k   EY dot V1_k   EZ dot V1_k
    EX dot V2_k   EY dot V2_k   EZ dot V2_k
    EX dot Vn_k   EY dot Vn_k   EZ dot Vn_k ]
```

The element-level local-to-global transformation is block diagonal:

```text
u_local_e = T_e u_global_e
K_global_e = T_e^T K_local_e T_e
f_global_e = T_e^T f_local_e
```

Each node block maps `[UX, UY, UZ, RX, RY, RZ]` to `[u_x, u_y, u_z, alpha, beta, gamma]`.

## Covariant Strain Definition
Use convected covariant components for the MITC4 strain construction:

```text
g_i = dX / dr_i, where r_i = xi, eta, zeta

eps_ij =
  0.5 * (du/dr_i dot g_j + du/dr_j dot g_i)
```

For the Phase 1 small-strain implementation:

- Membrane and bending terms are derived directly from the degenerated-continuum displacement field.
- `eps_33 = 0` is enforced by the shell assumption.
- The transverse shear rows `eps_13` and `eps_23` are not taken directly at the Gauss point. They are replaced by the MITC assumed transverse shear interpolation below.

Use the internal strain vector order:

```text
[ eps_11, eps_22, eps_33, gamma_23, gamma_13, gamma_12 ]^T
```

where engineering shear strains are:

```text
gamma_ij = 2 * eps_ij
```

If an implementation uses another row order internally, it must include an explicit permutation test against this documented order.

## MITC Transverse Shear Interpolation
The key MITC4 assumption is that the transverse shear tensor components are interpolated from midside tying points, instead of being evaluated directly at every Gauss point from the displacement field.

The original paper form is:

```text
eps_13 = 0.5 * (1 + r2) * eps_13^A + 0.5 * (1 - r2) * eps_13^C
eps_23 = 0.5 * (1 + r1) * eps_23^D + 0.5 * (1 - r1) * eps_23^B
```

With FESA's Abaqus-style natural coordinate convention and tying-point labels, use:

```text
eps_13_MITC(xi, eta) =
  0.5 * (1 - eta) * eps_13^A
  + 0.5 * (1 + eta) * eps_13^C

eps_23_MITC(xi, eta) =
  0.5 * (1 - xi) * eps_23^B
  + 0.5 * (1 + xi) * eps_23^D
```

The tying values are the direct covariant shear strains evaluated at midside points:

```text
eps_13^A = eps_13_DIRECT( 0, -1)
eps_13^C = eps_13_DIRECT( 0, +1)
eps_23^B = eps_23_DIRECT(-1,  0)
eps_23^D = eps_23_DIRECT(+1,  0)
```

Implementation rule:

- Build the direct covariant strain-displacement rows from the geometry and displacement interpolation.
- Evaluate only the needed direct transverse shear rows at `A`, `B`, `C`, and `D`.
- Interpolate those rows to the integration point using the formulas above.
- Do not compute MITC shear from local Cartesian `gamma_xz`/`gamma_yz` first; the tying is on convected tensor components.

For sign and regression checks, define:

```text
q_k = -V2_k * alpha_k + V1_k * beta_k
```

The paper edge formulas can be used as a diagnostic check for the direct tying rows:

```text
eps_13^A ~ 1/8 * [
  (u1 - u2) dot 0.5 * (t1 Vn1 + t2 Vn2)
  + (X1 - X2) dot 0.5 * (t1 q1 + t2 q2)
]

eps_13^C ~ 1/8 * [
  (u4 - u3) dot 0.5 * (t4 Vn4 + t3 Vn3)
  + (X4 - X3) dot 0.5 * (t4 q4 + t3 q3)
]

eps_23^B ~ 1/8 * [
  (u1 - u4) dot 0.5 * (t1 Vn1 + t4 Vn4)
  + (X1 - X4) dot 0.5 * (t1 q1 + t4 q4)
]

eps_23^D ~ 1/8 * [
  (u2 - u3) dot 0.5 * (t2 Vn2 + t3 Vn3)
  + (X2 - X3) dot 0.5 * (t2 q2 + t3 q3)
]
```

These compact formulas are not a license to skip the general direct-row evaluation. They are a sign/orientation cross-check derived from the same midside strain definitions.

## Local Cartesian Frame and Material Law
The constitutive law is defined in a local orthonormal Cartesian frame at each integration point with `sigma_33 = 0`.

Build the integration-point local basis from the covariant basis:

```text
e3_hat = normalize(g_3)
e1_hat = normalize(g_2 x e3_hat)
e2_hat = e3_hat x e1_hat
```

If `g_2 x e3_hat` is near zero, use a deterministic fallback that still creates a right-handed orthonormal basis and emits a diagnostic if no valid basis exists.

For isotropic linear elasticity, use the local plane-stress shell material matrix with transverse shear correction `kappa = 5/6`:

```text
D_hat = E / (1 - nu^2) *
  [ 1   nu  0   0                    0                    0
    nu  1   0   0                    0                    0
    0   0   0   0                    0                    0
    0   0   0   kappa*(1-nu)/2       0                    0
    0   0   0   0                    kappa*(1-nu)/2       0
    0   0   0   0                    0                    (1-nu)/2 ]
```

This matrix uses the documented strain order:

```text
[ eps_11, eps_22, eps_33, gamma_23, gamma_13, gamma_12 ]^T
```

For a general non-flat element, transform the local Cartesian material relation to the convected coordinate strain vector:

```text
E_hat = T_xi_to_x E_convected
D_convected = T_xi_to_x^T D_hat T_xi_to_x
```

The transformation matrix must be constructed from direction cosines between the local orthonormal basis and contravariant/covariant convected bases. For Phase 1, a flat-element implementation may compute the B-matrix directly in the local Cartesian frame only if tests prove it is equivalent to the convected formulation for the accepted reference cases.

## Element Stiffness and Numerical Integration
The baseline stiffness is:

```text
K_e = integral_V B^T D B dV
```

Phase 1 baseline integration:

- Use `2 x 2 x 2` Gauss integration in `(xi, eta, zeta)` for the first rebuilt MITC4 implementation.
- This follows the thesis derivation's `2 x 2 x 2` natural-coordinate integration and the Dvorkin-Bathe elastic recommendation of `2 x 2` in the midsurface with two through-thickness points.
- Do not introduce reduced integration or hourglass control.
- Analytic through-thickness integration may be introduced later only after equivalence tests preserve the MITC shear, drilling, and reference displacement behavior.

At every integration point:

1. Compute `g_i`, the metric/Jacobian, local Cartesian basis, and material transform.
2. Build direct membrane/bending strain rows.
3. Replace transverse shear rows with MITC interpolated rows.
4. Accumulate `B^T D B |J| w_xi w_eta w_zeta`.

Detect near-zero Jacobian, invalid basis, missing property/material, or non-positive thickness as invalid/singular element diagnostics. Do not report these as mesh quality diagnostics in Phase 1.

## Drilling DOF Stabilization
The continuum MITC4 formulation has five physical DOFs per node. FESA keeps six global DOFs, so the local drilling rotation `gamma` has no physical strain contribution and must be weakly stabilized.

Baseline decision:

```text
drilling_stiffness_scale = 1.0e-3
k_ref = min_positive_diagonal(K_local_without_drilling)
k_drill = drilling_stiffness_scale * k_ref
```

Rules:

- `K_local_without_drilling` is the element local stiffness after MITC integration and before artificial drilling stabilization.
- Use the minimum positive diagonal entry over the physical local DOFs. Ignore zero, negative, NaN, or non-finite entries.
- Add `k_drill` to each nodal local `gamma` diagonal before transforming to global coordinates.
- If no positive reference diagonal exists, emit a singular/invalid element diagnostic.
- The scale must be configurable and written to result metadata or analysis diagnostics.
- Reference comparisons must include a drilling sensitivity check if displacement results change materially with the scale.

This replaces the earlier arbitrary `1.0e-6 * E * thickness` rule. The `1.0e-3` scale follows the local thesis derivation's six-DOF stabilization note and is still an artificial numerical parameter, not a physical stiffness.

## Internal Force and Reaction Recovery
For Phase 1 linear static analysis:

```text
f_int_e = K_e u_e
R_full = K_full U_full - F_full
```

Rules:

- Reactions must be recovered from the full global vector, not only from the reduced free-DOF system.
- `RF` output follows the Abaqus-style component convention documented in `docs/NUMERICAL_CONVENTIONS.md`.
- Element stress/strain/resultant outputs are not mandatory until displacement and reaction reference behavior is stable.

## Required Tests Before Reimplementation
The rebuilt MITC4 element must be protected by tests before implementation code is accepted.

Element-level tests:

- Shape functions sum to one and have correct derivatives at center, corners, and tying points.
- Center normal, local director axes, and integration local Cartesian axes are orthonormal and right-handed.
- The local-to-global rotation transform maps global rotations to `alpha`, `beta`, `gamma` as documented.
- Direct tying rows at `A`, `B`, `C`, `D` agree with finite-difference strain perturbations.
- MITC shear rows at Gauss points are interpolated from tying rows, not evaluated directly at the Gauss points.
- Stiffness is symmetric for linear elastic Phase 1.
- Rigid-body translations and rotations produce near-zero physical strain energy, with only documented drilling stabilization effects.
- The element has no spurious zero-energy modes beyond expected rigid-body behavior after drilling stabilization is accounted for.

Patch and benchmark tests:

- Constant membrane states.
- Pure bending.
- Pure shear.
- Pure twist.
- Thin cantilever or plate strip to expose shear locking.
- Scordelis-Lo roof after the basic flat tests pass.
- Pinched cylinder or twisted beam as later Phase 1/Phase 2 benchmark expansion.

Reference regression tests:

- Compare FESA `U` against stored Abaqus displacement CSV files using absolute and relative tolerances.
- Keep `RF` verified by full-vector equilibrium until matching Abaqus reaction CSV files are provided or explicitly deferred.

## Current Reference Artifacts
Stored artifacts currently known:

- `references/quad_01.inp`
- `references/quad_01_displacements.csv`
- `references/quad_02.inp`
- `references/quad_02_displacements.csv`

Compatibility notes:

- `quad_01.inp` contains `S4R`, `Part/Assembly/Instance`, `*Density`, and `NLGEOM=YES`; it remains future compatibility provenance, not a Phase 1 parser acceptance case.
- `quad_02.inp` contains `TYPE=S4`, which is the correct element target for Phase 1, but it also uses `Part/Assembly/Instance`. The next Phase 1 planning pass must either normalize this input into the Phase 1 flat keyword subset or explicitly add a parser sprint for this Abaqus/CAE structure. Do not silently expand parser support while implementing the element.

## Implementation Checklist
The next MITC4 implementation pass should proceed in this order:

1. Lock the shape-function, natural-coordinate, node-order, and tying-point tests.
2. Implement geometry/director/local-axis construction with diagnostics.
3. Implement the five-DOF degenerated-continuum displacement interpolation.
4. Implement global six-DOF to local `[alpha, beta, gamma]` transformation.
5. Implement direct covariant strain rows and finite-difference tests.
6. Replace transverse shear rows with MITC tying interpolation.
7. Implement material transform and `2 x 2 x 2` stiffness integration.
8. Add drilling stabilization in local coordinates.
9. Transform element stiffness and internal force to global coordinates.
10. Run element patch tests before assembly/reference regression work.

## Future Extensions
Geometric nonlinearity:

- Update current geometry and nodal directors in `AnalysisState`.
- Add finite-rotation updates for `V1`, `V2`, and `Vn`.
- Add tangent stiffness, geometric stiffness, and Newton-Raphson convergence checks.

Thermal-stress coupling:

- Add temperature field state and thermal strain contribution.
- Add material expansion data and thermal output fields.

S4R:

- Add only after a separate ADR and formulation update.
- Requires reduced integration, hourglass control, and separate reference artifacts.

Stress/resultant output:

- Define section point locations, component labels, local coordinate convention, and Abaqus CSV comparison format before making `S`, `E`, or `SF` mandatory.

## Remaining Open Decisions
The MITC4 stiffness formulation is now defined for the Phase 1 rebuild. Remaining decisions are outside the first stiffness pass:

- Exact stress/strain/resultant recovery locations and component labels.
- Whether to parse Abaqus `*Orientation` or nodal shell normals.
- Whether `quad_02.inp` should be normalized into the existing Phase 1 parser subset or used to justify a dedicated `Part/Assembly/Instance` parser sprint.
- Reference tolerances for `quad_02` after the input compatibility path is chosen.