baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

docs: reset mitc4 formulation baseline

Browse Source
This commit is contained in:
NINI
2026-05-04 12:00:16 +09:00
parent c5be1f988c
commit a385235e14
117 changed files with 5924 additions and 213 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/MITC4_FORMULATION.md
+375 -191
View File
@@ -1,257 +1,441 @@
# MITC4 Formulation
## Purpose
This document defines the baseline MITC4 formulation target for FESA Phase 1.
This document is the formulation contract for FESA's Phase 1 four-node MITC4 shell element.
It is intentionally a formulation contract, not implementation code. Exact formulas should be added and reviewed before coding the MITC4 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.
## 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
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.
## Phase 1 Target
Phase 1 implements a clear MITC4 baseline formulation and passes reference benchmarks before performance optimization.
## Local Source Basis
The formulation below is based on the local papers in `docs/Paper/`:
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.
- `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.
Non-scope:
- S4R reduced-integration behavior.
- Hourglass control.
- Composite sections.
- Material nonlinearity.
- Geometric nonlinearity.
- Pressure loads.
- Thermal-stress coupling.
- Mesh quality diagnostics.
## Phase 1 Scope
Phase 1 implements:
## 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.
- 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`.
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:
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
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)
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
```
- 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:
Shape functions:
```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))
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)
```
- 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:
Midside MITC tying points in the FESA convention:
```text
drilling_stiffness_scale = 1.0e-6
A = ( 0, -1) edge 1-2
B = (-1, 0) edge 1-4
C = ( 0, +1) edge 4-3
D = (+1, 0) edge 2-3
```
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.
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.
## Nodal DOFs
Each node has:
## Degenerated-Continuum Geometry
The shell is represented as a 3D continuum degenerated through the thickness:
```text
UX, UY, UZ, RX, RY, RZ
X(xi, eta, zeta) =
sum_k N_k(xi, eta) X_k
+ zeta / 2 * sum_k t_k N_k(xi, eta) Vn_k
```
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.
where:
## 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.
- `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.
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.
Phase 1 director policy:
## 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:
- If no parsed nodal directors exist, use the element-center midsurface normal for all four nodal directors:
```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)
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
```
where `r, s` are natural coordinates in `[-1, 1]`.
- 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.
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.
The nodal local director axes are:
## Strain Treatment
The baseline element separates:
- membrane strain terms.
- bending curvature terms.
- transverse shear strain terms.
- artificial drilling stabilization.
```text
V1_k = normalize(EY x Vn_k)
V2_k = Vn_k x V1_k
```
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.
If `EY x Vn_k` is near zero, use a deterministic fallback axis before normalization:
The Phase 1 transverse shear tying formula is the midside interpolation stated in `Phase 1 Closed Baseline Decisions`.
```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
```
## 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.
The implementation must keep `V1_k`, `V2_k`, and `Vn_k` orthonormal and right-handed.
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.
## 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
Decision:
- Phase 1 uses small artificial drilling stiffness.
- Default scale: `drilling_stiffness_scale = 1.0e-6`.
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.
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:
Baseline decision:
```text
k_drill = drilling_stiffness_scale * representative_element_stiffness
representative_element_stiffness = E * thickness
drilling_stiffness_scale = 1.0e-3
k_ref = min_positive_diagonal(K_local_without_drilling)
k_drill = drilling_stiffness_scale * k_ref
```
The representative stiffness may be refined after the first Abaqus S4 reference cases are available.
Rules:
## 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.
- `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.
Future output:
- local shell stresses.
- local shell strains.
- section forces and moments.
- integration point and section point data.
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.
## 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.
## Internal Force and Reaction Recovery
For Phase 1 linear static analysis:
## 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.
```text
f_int_e = K_e u_e
R_full = K_full U_full - F_full
```
Stress/strain recovery locations for `S`, `E`, and `SF` remain future decisions and must not be implemented as mandatory Phase 1 outputs.
Rules:
## 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.
- 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.
Distorted mesh tests should be added after the baseline passes, but Phase 1 does not implement general mesh quality diagnostics.
## 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:
Current stored reference artifact:
- `references/quad_01.inp`
- `references/quad_01_displacements.csv`
- `references/quad_02.inp`
- `references/quad_02_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.
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:
- Add updated geometry, current frame handling, tangent stiffness, and Newton-Raphson integration.
- Preserve `AnalysisState` element/internal state extension points.
- 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.
- Add thermal strain contribution.
- Add material expansion data.
- Add result fields for temperature and thermal strain/stress.
- 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 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.
- 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.