FESADev/docs/MITC4_FORMULATION.md

# MITC4 Formulation

## Purpose
This document defines the baseline MITC4 formulation target for FESA Phase 1.

It is intentionally a formulation contract, not implementation code. Exact formulas should be added and reviewed before coding the MITC4 element.

## Source Basis
- Dvorkin and Bathe's four-node shell element paper presents a continuum-mechanics-based, non-flat, general quadrilateral shell element for thin and thick shells and nonlinear analysis: https://web.mit.edu/kjb/www/Publications_Prior_to_1998/A_Continuum_Mechanics_Based_Four-Node_Shell_Element_for_General_Nonlinear_Analysis.pdf
- The paper identifies transverse shear locking as a key problem in simple 4-node shell interpolation and motivates modified transverse shear treatment: https://web.mit.edu/kjb/www/Publications_Prior_to_1998/A_Continuum_Mechanics_Based_Four-Node_Shell_Element_for_General_Nonlinear_Analysis.pdf
- OpenSees describes `ShellMITC4` as a bilinear isoparametric shell element with modified shear interpolation, four counter-clockwise nodes, and six DOFs per node: https://opensees.berkeley.edu/wiki/index.php/Shell_Element
- The MITC benchmark paper states that the MITC method is used to remedy shell locking and that the standard MITC4 employs MITC treatment for transverse shear strains; it also notes that Abaqus S4 uses Dvorkin-Bathe transverse shear interpolation: https://web.mit.edu/kjb/www/Principal_Publications/Performance_of_the_MITC3%2B_and_MITC4%2B_shell_elements_in_widely_used_benchmark_problems.pdf
- Abaqus finite-strain shell theory documentation provides useful comparison context for S4/S4R geometry, interpolation, orientation update, and transverse shear treatment, but FESA Phase 1 is linear static: https://abaqus-docs.mit.edu/2017/English/SIMACAETHERefMap/simathe-c-finitestrainshells.htm

## Phase 1 Target
Phase 1 implements a clear MITC4 baseline formulation and passes reference benchmarks before performance optimization.

Scope:
- 4-node quadrilateral shell.
- Linear static analysis.
- Linear isotropic elastic material.
- Homogeneous shell section.
- 6 DOFs per node.
- Small-strain formulation for Phase 1.
- Transverse shear interpolation based on MITC4 assumptions.
- Abaqus-compatible result signs.

Non-scope:
- S4R reduced-integration behavior.
- Hourglass control.
- Composite sections.
- Material nonlinearity.
- Geometric nonlinearity.
- Pressure loads.
- Thermal-stress coupling.
- Mesh quality diagnostics.

## Phase 1 Closed Baseline Decisions
The following decisions close the Phase 1 implementation gate for the first baseline C++ implementation. They are intentionally conservative and may be revised by ADR after Abaqus S4 reference cases are added.

Decisions:
- Mandatory Phase 1 result outputs are limited to nodal `U` and `RF`. Element `S`, `E`, and `SF` remain future outputs.
- In-plane membrane and bending terms use standard bilinear Reissner-Mindlin interpolation with 2x2 Gauss integration.
- Transverse shear uses MITC4 mid-side tying interpolation:

```text
gamma_xz_mitc(r, s) =
  0.5 * (1 - s) * gamma_xz(0, -1)
  + 0.5 * (1 + s) * gamma_xz(0, +1)

gamma_yz_mitc(r, s) =
  0.5 * (1 - r) * gamma_yz(-1, 0)
  + 0.5 * (1 + r) * gamma_yz(+1, 0)
```

- The four tying points are the midside natural-coordinate points `(0,-1)`, `(0,+1)`, `(-1,0)`, and `(+1,0)`.
- Local basis construction uses averaged opposite edges:

```text
v1 = 0.5 * ((x2 - x1) + (x3 - x4))
v2 = 0.5 * ((x4 - x1) + (x3 - x2))
e1 = normalize(v1)
e2_raw = v2 - dot(v2, e1) * e1
e2_pre = normalize(e2_raw)
e3 = normalize(cross(e1, e2_pre))
e2 = normalize(cross(e3, e1))
```

- For mildly warped quadrilaterals this is an averaged midsurface basis. Severe warpage is not diagnosed as mesh quality in Phase 1, but near-zero vectors or Jacobians are invalid/singular element diagnostics.
- Integration point order for stiffness tests is `(-g,-g)`, `(g,-g)`, `(g,g)`, `(-g,g)` conceptually, with `g = 1 / sqrt(3)`. Implementation may loop over the same set in row-major order as long as assembled stiffness is invariant.
- The default drilling stiffness scale is:

```text
drilling_stiffness_scale = 1.0e-6
```

The current baseline applies this scale to a representative `E * thickness` drilling stabilization term. The value is reported through the element options path and must be revisited after Abaqus S4 displacement references are available.

## Nodal DOFs
Each node has:

```text
UX, UY, UZ, RX, RY, RZ
```

Rules:
- Translational DOFs are global translations.
- Rotational DOFs are rotations about global or transformed axes following Abaqus component convention.
- `RZ` is retained as a drilling DOF.
- Drilling stiffness is artificial in Phase 1 and must be parameterized.

## Element Input Contract
`MITC4Element` requires:
- four node ids in Abaqus S4 order.
- four node coordinates.
- shell thickness.
- linear elastic material constants `E` and `nu`.
- drilling stiffness parameter.
- element id and property id for diagnostics and output.

Node ordering:
- Input node order follows Abaqus S4 convention.
- Positive normal follows the right-hand rule around the nodes.
- FESA maps Abaqus `TYPE=S4` to `MITC4`.
- Abaqus `TYPE=S4R` is not supported in Phase 1.

## Coordinate Frames
The Phase 1 local basis construction is defined by the averaged-edge algorithm in `Phase 1 Closed Baseline Decisions`.

Minimum requirements:
- Define a local shell normal from the quadrilateral geometry.
- Define local in-plane axes `e1` and `e2` so that `e1`, `e2`, and normal form a right-handed basis.
- Preserve Abaqus-compatible output signs.
- Document behavior for non-planar quadrilaterals.
- Use the same convention consistently for stiffness, stress/strain recovery, and result output.

Phase 1 convention:
- Use the element midsurface geometry to compute an average basis from opposite edges.
- Use the positive shell normal implied by the input node ordering.
- Use a right-handed local basis consistently for stiffness, reaction recovery, and future result output.

## Shape Functions
Baseline quadrilateral bilinear interpolation:

```text
N1 = 0.25 * (1 - r) * (1 - s)
N2 = 0.25 * (1 + r) * (1 - s)
N3 = 0.25 * (1 + r) * (1 + s)
N4 = 0.25 * (1 - r) * (1 + s)
```

where `r, s` are natural coordinates in `[-1, 1]`.

Implementation requirements:
- Compute shape function derivatives with respect to natural coordinates.
- Build the surface Jacobian and local derivatives.
- Detect invalid or near-zero Jacobian as a singular/invalid element diagnostic, not as a mesh quality metric.

## Strain Treatment
The baseline element separates:
- membrane strain terms.
- bending curvature terms.
- transverse shear strain terms.
- artificial drilling stabilization.

MITC4 requirement:
- Use standard displacement interpolation for membrane and bending terms in Phase 1.
- Use MITC transverse shear interpolation to alleviate shear locking.
- Do not replace MITC4 with plain full-integration Reissner-Mindlin Q4.

The Phase 1 transverse shear tying formula is the midside interpolation stated in `Phase 1 Closed Baseline Decisions`.

## Numerical Integration
Phase 1 baseline:
- In-plane integration: 2x2 Gauss for membrane, bending, MITC shear, and drilling stabilization.
- Thickness integration: homogeneous linear elastic shell section is integrated analytically using `t`, `t^3 / 12`, and shear correction `5 / 6`.
- Benchmark literature commonly reports 2x2 in-plane Gauss integration for S4/MITC4-style 4-node elements and 2-point thickness integration in comparative shell studies.

Rules:
- Do not introduce reduced integration or hourglass control for S4R behavior in Phase 1.
- Do not optimize integration before reference benchmarks pass.
- Integration point ordering for output must be documented before stress/strain reference comparisons.

## Drilling DOF Stabilization
Decision:
- Phase 1 uses small artificial drilling stiffness.
- Default scale: `drilling_stiffness_scale = 1.0e-6`.

Requirements:
- Expose a parameter such as `drilling_stiffness_scale`.
- Provide a deterministic default.
- Make the default small enough not to dominate physical response.
- Include benchmark sensitivity checks if reference results are sensitive to the value.
- Report the value in result metadata.

Baseline rule:

```text
k_drill = drilling_stiffness_scale * representative_element_stiffness
representative_element_stiffness = E * thickness
```

The representative stiffness may be refined after the first Abaqus S4 reference cases are available.

## Element Outputs
Phase 1 minimum:
- element stiffness matrix.
- element equivalent nodal internal force for full-vector residual/reaction recovery.
- optional stress/strain output after displacement benchmarks are stable.

Future output:
- local shell stresses.
- local shell strains.
- section forces and moments.
- integration point and section point data.

## Required Element-Level Tests
Before integration with the global solver:
- shape functions sum to one.
- derivatives satisfy expected bilinear identities.
- element stiffness dimensions are `24 x 24`.
- stiffness is symmetric for linear elastic Phase 1.
- rigid body translations produce near-zero internal strain energy.
- rigid body rotations do not create physical membrane/bending stiffness beyond documented drilling effects.
- constant membrane patch behavior.
- bending-dominated sanity case.
- drilling stiffness sensitivity check.

## Pre-Implementation Gate
Phase 1 baseline implementation may proceed because these items are now recorded above:
- transverse shear tying point equations.
- local shell basis algorithm for flat and mildly warped quadrilaterals.
- default drilling stiffness scale and parameter name.
- integration point set and ordering convention.
- Phase 1 mandatory output scope: `U` and `RF` only.

Stress/strain recovery locations for `S`, `E`, and `SF` remain future decisions and must not be implemented as mandatory Phase 1 outputs.

## Reference Benchmarks
MITC4 baseline acceptance should include:
- single-element membrane test.
- single-element bending test.
- cantilever shell strip.
- simply supported square plate.
- Scordelis-Lo roof.
- pinched cylinder.
- twisted beam.

Distorted mesh tests should be added after the baseline passes, but Phase 1 does not implement general mesh quality diagnostics.

Current stored reference artifact:
- `references/quad_01.inp`
- `references/quad_01_displacements.csv`

Compatibility note:
- `quad_01.inp` is an Abaqus-generated S4R/NLGEOM reference input and is not a Phase 1 MITC4 parser acceptance case.
- It should be used as an external reference artifact and future compatibility target unless a normalized Phase 1 `S4` input is added.
- The displacement CSV can still exercise `U` field comparison infrastructure once node ids and component mapping are supported.

## Future Extensions
Geometric nonlinearity:
- Add updated geometry, current frame handling, tangent stiffness, and Newton-Raphson integration.
- Preserve `AnalysisState` element/internal state extension points.

Thermal-stress coupling:
- Add temperature field state.
- Add thermal strain contribution.
- Add material expansion data.
- Add result fields for temperature and thermal strain/stress.

S4R:
- Add only after a separate ADR and formulation document update.
- Requires reduced integration and hourglass control decisions.

## Open Decisions Before Coding
- Exact stress/strain recovery locations.
- Whether local coordinate transforms from Abaqus input are deferred or rejected.