# OpenSees ShellMITC4 quad results vs CalculiX S3 **Date:** 2026-05-04 **Author:** Claude (Phase 2.3 subagent) **Scope:** Add a true Mindlin-Reissner shell route (OpenSees ShellMITC4 on quads) to the dome analysis stack and compare against the trusted CalculiX S3 reference. Companion to `reports/2026-05-04--shell-solver-disagreement.md`. --- ## TL;DR ShellMITC4 (Mindlin-Reissner quad) does **NOT** agree with CCX S3 (Mindlin-Reissner triangle) on the curved dome. On the same load case (`wind_cc_peak`, baseline) MITC4 gives `u_max = 29.83 mm` (4 sub-quads per panel) and `u_max = 36.64 mm` (64 sub-quads per panel), versus **CCX S3's 12.78 mm**. That is a 2.3-2.9× over-prediction — about the same band as the in-house DKT result. **Mesh refinement makes it worse, not better.** The flat-plate Timoshenko smoke test passes: MITC4 on a 10×10 quad clamped square gives `ratio_FEA/Kirchhoff = 1.108`, in the expected Mindlin-Reissner shear-deformation band [1.00, 1.20] for `a/t = 13.1` (DKT ratio's 1.024 on the same mesh — DKT misses the Mindlin shear contribution by construction). So the formulation itself is wired correctly on flat geometry. **The disagreement is specific to the curved dome**, exactly as it was for the two DKT solvers. This is the **second** solver-disagreement finding. Phase 2 still has one trusted Mindlin-Reissner shell route (CCX S3); the new MITC4 path behaves more like the DKT family than like CCX. Hypotheses on the root cause are listed in §5; further investigation is required before MITC4 can be used as a Phase-3/4 shell route. --- ## 1. Comparison table — wind_cc_peak baseline Site: baseline CONUS. Load case: `wind_cc_peak` (p = −3 831 Pa, outward suction). Same mesh strategy (same midsurface corner network, `dedup_tol_mm = 300`) is used for the in-house DKT, the OpenSees ShellDKGT, and the new OpenSees ShellMITC4 runs. CCX S3 uses the same 8 960-triangle midsurface mesh. | Solver | u_max (mm) | max σ_b (MPa) | Formulation | Mesh | |---|---|---|---|---| | **CalculiX S3** | **12.78** | n/a (CCX OUTPUT=2D, membrane only) | Mindlin-Reissner triangular shell | 8 960 tri, 4 639 nodes | | OpenSees ShellMITC4 (4 sub-quads/panel) | 29.83 | 0.446 | Mindlin-Reissner quad shell | 280 quad, 319 nodes | | OpenSees ShellMITC4 (64 sub-quads/panel)* | 36.64 | (not aggregated)| Mindlin-Reissner quad shell | 4 480 quad, 4 639 nodes | | In-house CST + DKT | 29.47 | 0.645 | Discrete Kirchhoff Triangle | 8 960 tri, 4 639 nodes | | OpenSees ShellDKGT | 52.42 | 1.745 | Discrete Kirchhoff Triangle | 8 960 tri, 4 639 nodes | \* mesh-refinement experiment, not the production sweep — see §3. CCX σ_b is left blank because the existing CCX shell run was configured with `OUTPUT=2D` which writes only membrane (in-plane) stress to the .frd, not bending moments. (`reports/ccx_shell/baseline_wind_cc_peak.frd` contains 4 639 nodal stress tensors; running it with `OUTPUT=3D` to recover bending stresses is a follow-up.) --- ## 2. Smoke-test results `tools/opensees_shell_mitc_smoke.py` (10×10 quads, clamped square, PU-foam material, q = 1000 Pa): ``` w_centre (FEA) : 0.4866 mm w_centre (theory) : 0.4392 mm (Kirchhoff, Timoshenko Table 35) ratio FEA/theory : 1.1079 ``` The 11% over-Kirchhoff is *expected* for a Mindlin-Reissner element on this slenderness (`a/t = 13.1`). It is the transverse shear contribution to centre deflection — DKT (Discrete Kirchhoff) would not produce this offset because it omits transverse shear by construction. CCX S3 will produce the same Mindlin-vs-Kirchhoff offset on the same problem. **Verdict on the smoke test: PASS.** The ShellMITC4 wrapper is wired correctly on flat geometry. --- ## 3. Production sweep — `tools/solve_shell_quad_opensees.py baseline` Run output: `reports/full-analysis-shell-quad-opensees-baseline.txt`. Per-panel CSV: `reports/opensees_shell_quad_baseline_forces.csv` (70 panels × 3 cases = 210 rows). | Load case | p (Pa) | u_max (mm) | σ_b_max (MPa) | σ_m_max (MPa) | VM_max (MPa) | D/C bending | |---|---|---|---|---|---|---| | snow | +1 005 | 8.55 | 0.113 | 0.235 | 0.315 | 0.13 | | wind_mwfrs | −1 419 | 11.14 | 0.167 | 0.102 | 0.227 | 0.19 | | wind_cc_peak | −3 831 | 29.83 | 0.446 | 0.400 | 0.766 | 0.51 | For comparison, the existing OpenSees DKT triangle run on the same site gives wind_cc_peak `u_max = 52.42 mm` and σ_b_max = 1.745 MPa (`reports/full-analysis-shell-opensees-baseline.txt`). MITC4's quad result is **smaller than DKT triangle**, but still **2.33× larger than CCX S3**. Top-10 hottest panels are `64, 62, 10, 2, 63, 59, 58, 61, 60, 48` — same regions both solvers identify as hottest, just at different absolute magnitudes. ### 3a. Mesh-refinement check Refining each panel into N×N sub-quads (rather than the 2×2 = 4 sub-quads of the production sweep) under the same load: | n_div | quads | nodes | u_max (mm) | ratio to CCX | |---|---|---|---|---| | 2 | 280 | 319 | 29.83 | 2.33 | | 4 | 1 120 | 1 199 | 31.22 | 2.44 | | 6 | 2 520 | 2 639 | 35.00 | 2.74 | | 8 | 4 480 | 4 639 | 36.64 | 2.87 | The n_div=8 quad mesh has the **same 4 639 nodes** as the CCX S3 triangle mesh (because the corner network and edge subdivision are both driven by `dedup_tol_mm=300` and `target_size_mm=200`). Same nodes, same BCs, same load lumping convention — yet CCX gives 12.78 mm and OpenSees MITC4 gives 36.64 mm. Refinement makes the disagreement **larger**, not smaller, so it is **not** a coarse-mesh shear-locking artefact (which would converge *down* to the true answer with refinement). The behaviour is qualitatively different from the standard MITC4 known issue of membrane locking on coarse curved meshes (which would *under*-predict deflection). --- ## 4. Comparison versus the prior 3-solver picture Prior (`reports/2026-05-04--shell-solver-disagreement.md`): | Solver | u_max (mm) | Ratio to CCX S3 | |---|---|---| | CalculiX S3 | 12.78 | 1.00 | | In-house CST + DKT | 29.47 | 2.31 | | OpenSees ShellDKGT | 52.42 | 4.10 | Now (with new MITC4 row): | Solver | u_max (mm) | Ratio to CCX S3 | Formulation | |---|---|---|---| | CalculiX S3 | 12.78 | 1.00 | Mindlin-Reissner tri | | **OpenSees ShellMITC4 (this work)** | **29.83** | **2.33** | **Mindlin-Reissner quad** | | In-house CST + DKT | 29.47 | 2.31 | Kirchhoff tri | | OpenSees ShellDKGT | 52.42 | 4.10 | Kirchhoff tri | MITC4 stacks essentially **on top of the in-house DKT** in this band, not on top of CCX. The "physics-class" hypothesis from the disagreement report (DKT-vs-Mindlin explains the 2-4× spread) is now falsified: a true Mindlin-Reissner shell on the same dome still gives ~2.3× over CCX. So the gap **cannot** be attributed solely to DKT's missing shear and missing membrane-bending coupling. --- ## 5. Hypotheses (ordered by what to investigate first) These are **not yet diagnosed** — listing in priority order so a follow-up investigation has a starting plan. 1. **CCX S3 stiffness recipe.** CalculiX's S3 element is a 3-node triangle but its formulation is more elaborate than a stock Mindlin-Reissner tri — it uses an MITC-style shear interpolation and may include element-level drilling stabilization. If CCX S3 is *stiffer* than the textbook Mindlin-Reissner shell on faceted curved geometry, then the apparent "MITC4 over-predicts" gap may actually be "CCX S3 under-predicts because of an extra stabilization term we are not matching." Step: cross-check CCX against an independent S3 result (Code_Aster shells, Abaqus S3R). 2. **Quad outward-normal alignment.** `build_quad_midsurface_mesh` flips 100/280 quad windings to make the right-hand-rule normal match the panel's `FittedPanel.normal`. If the panel-level normal convention is itself slightly off (e.g. averaged-normal vs PCA-best-fit-normal mismatch), the pressure load lumping per quad could be biased. Verify by recomputing per-quad outward normal from the 4 corner positions only (no reference to panel normal), re-lumping, and comparing. 3. **Per-panel non-planarity.** Each rhombic panel's 4 cluster centroids are intentionally NOT perfectly coplanar (averaged normals across joints introduce small out-of-plane offsets). The triangle midsurface mesh handled this by splitting each panel into 2 triangles along the rhombus diagonal, which are guaranteed planar. The quad mesh splits each panel into 4 sub-quads which are *bilinearly interpolated* on the 4 corner positions — they are NOT guaranteed planar, and ShellMITC4 may not handle warped quads as gracefully as CCX S3 handles warped triangles. Step: measure per-quad warpage (max corner-to-best-fit-plane distance in mm) and correlate with disagreement. 4. **Shear-strain interpolation on a coarse curved mesh.** MITC4's "mixed interpolation of tensorial components" for transverse shear was developed and validated on smooth curved shells with smooth meshes. A faceted dome with 4 sub-quads per facet has a very abrupt curvature jump at every panel-to-panel joint. The MITC shear interpolation may over-respond to those joint-line discontinuities. Step: replace the faceted dome with a smooth spherical dome of similar diameter, mesh with quads of similar density, run MITC4 — if it agrees with CCX on the smooth dome but not the faceted one, this hypothesis stands. 5. **Drilling DOF / rotational coupling.** ShellMITC4 uses a penalty drilling stiffness; tweaking it might change behaviour on faceted curved shells. Mirror the in-house drilling sweep from `tools/investigate_drilling_sweep.py`. --- ## 6. Recommendation for Phase 3 / Phase 4 **Do NOT promote OpenSees ShellMITC4 to a trusted shell solver yet.** The 2.33× disagreement vs CCX S3 on the production load case is too large to ignore, and refinement makes it worse. Continue treating **CCX S3 as the primary trusted shell solver** for the dome (as already adopted in §6 of the disagreement report). The MITC4 wrapper is correctly wired (smoke test passes) and is a useful artifact: - It can serve as a future-investigation hook for the §5 hypotheses above. - It produces ranked hot-panel orderings consistent with both CCX and in-house DKT (panels 64, 62, 10, 2, 63, 59 always lead) — i.e. the *qualitative* stress map is right even if the absolute magnitudes are off, so it is still useful as a sanity check. - It can be re-enabled as a real third route once one of the §5 hypotheses is resolved. **Action items for Phase 3 (buckling) and Phase 4 (code-check):** - Phase 3: use CCX S3 (already validated) as the buckling-stiffness source. Treat in-house DKT and OpenSees DKT/MITC4 as conservative upper-bound sanity checks. - Phase 4: code-check using CCX S3 stresses. Document MITC4 as a known-disagreement that requires investigation before promotion. If the §5 investigation lands successfully, MITC4 can become the second trusted route in a future iteration. --- ## 7. Files Implementation (this Phase 2.3 work): - `src/zomestruct/fea/quad_mesh.py` — new - `src/zomestruct/fea/opensees_shell.py` — added `solve_static_shell_quad` - `tools/opensees_shell_mitc_smoke.py` — new flat-plate smoke - `tools/solve_shell_quad_opensees.py` — new production sweep Outputs: - `reports/opensees_shell_mitc_smoke.txt` — smoke pass - `reports/full-analysis-shell-quad-opensees-baseline.txt` — production run - `reports/opensees_shell_quad_baseline_forces.csv` — per-panel forces - `reports/opensees_shell_quad_baseline_*.vtu` — per-case VTUs (not committed) Reference results (already on disk): - `reports/ccx_shell/baseline_wind_cc_peak.frd` — CCX S3 - `reports/full-analysis-shell-opensees-baseline.txt` — OpenSees DKT - `reports/2026-05-04--shell-solver-disagreement.md` — prior 3-solver finding --- ## Update — smooth-cap test (2026-05-04 evening) A follow-up 4-way comparison on a **smooth spherical cap** of the same overall geometry as the production dome (R = 2950 mm, apex H = 3820 mm, θ_rim = 107.15°, t = 76.2 mm) under the same wind_cc_peak baseline load shows that **all four shell solvers agree within ±7%** on the smooth geometry: CCX S3 = 2.588 mm, MITC4 = 2.766 mm (ratio 1.069), in-house DKT = 2.581 mm (ratio 0.997), OpenSees DKGT = 2.628 mm (ratio 1.016). This **confirms the §5.4 hypothesis** (and falsifies §5.2 and most of §5.3) — the 2.33× CCX-vs-MITC4 disagreement on the production dome is caused by **faceting** (sharp-angle joints between flat panels), not by any intrinsic difference between CCX's and MITC4's Mindlin-Reissner formulations on smooth curved geometry. CCX S3 evidently has special handling for faceted curved shells that the other three solvers lack (or vice-versa — direction not yet pinned down). See [`2026-05-04--smooth-cap-4way.md`](2026-05-04--smooth-cap-4way.md) for the full table, mechanism discussion, and revised Phase 3 / Phase 4 recommendation (multi-solver D/C envelope).