# Shell-solver disagreement on the curved dome **Date:** 2026-05-04 **Author:** Claude (investigation subagent) **Scope:** Investigate why `OpenSees ShellDKGT` and the in-house `CST+DKT` shell solver disagree by ~2× on the curved dome under wind_cc_peak, yet agree to ~5% on a flat clamped plate. **Steps run:** A (re-confirmation, from existing artifacts) + B (per-node disagreement map) + C (drilling-stiffness sweep) + D (CCX S3 tiebreaker). Step E (mesh-refinement study) deferred. --- ## TL;DR The disagreement is **diffuse across the interior of the mesh** (median fractional disagreement on interior nodes ~28%, on joint nodes ~26%, on boundary nodes 0%). It is therefore **not** a boundary-condition or drilling-DOF artifact. A drilling-factor sweep on the in-house solver (0 → 1.0) moves `u_max` by less than ~22% — never approaching the OpenSees value of 52 mm (in-house at every drilling factor lies in the 24–30 mm band). Drilling stabilization can be ruled out as the root cause. The CalculiX S3 shell — independent third reference, same mesh — gives `u_max = 12.78 mm`. Under wind_cc_peak the three solvers stack like: | 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 | So the ~2× in-house ↔ OpenSees gap is part of a larger ~4× spread between **all three** shell solvers on the curved dome, even though all three pass the flat-plate Timoshenko reference within ~5%. The most likely root cause is the well-known **inadequacy of pure DKT (Discrete Kirchhoff Triangle) elements on faceted approximations of curved shells**: DKT is a Kirchhoff plate-bending element with no transverse shear and no membrane–bending coupling at the element level, so it sees inter-facet curvature only as discrete folds across element edges. The two DKT implementations differ in how they construct the local frame and rotate to global coords on each non-coplanar joint edge, which produces the 1.8× spread between in-house and OpenSees. CCX S3 includes shear coupling and membrane–bending interaction, so it recovers the curved-shell stiffness much better. **Recommendation:** treat CCX S3 as the trusted reference for the dome, treat both pure-DKT solvers (in-house, OpenSees) as conservative upper bounds, and **for Phase 3 buckling and Phase 4 code-check do NOT proceed with all three in parallel as Option C suggested** — instead adopt CCX S3 as the primary shell solver and use the OpenSees / in-house DKT runs only as an upper-bound sanity check. See § 6 below. --- ## 1. The disagreement (re-stated with current numbers) Site: baseline CONUS. Load case: `wind_cc_peak` (p = −3 831 Pa, outward suction). Shared inputs: - Mesh: `MidSurfaceMesh` from `build_midsurface_mesh(target_size_mm=200)` → 4 639 nodes, 8 960 triangles, 70 panels, 338 boundary nodes. - Material: PuFoam (E = 70.8 MPa, ν = 0.30, ρ = 240 kg/m³, t = 76.2 mm). - BC: encastre at every boundary node (all 6 DOFs). - Load: uniform per-triangle pressure with `area * p / 3` lumped to each corner translation DOF along the outward normal. | Solver | u_max [mm] | σ_b_max [MPa] | σ_m_max [MPa] | VM_max [MPa] | |---|---|---|---|---| | CalculiX S3 (`reports/ccx_shell/baseline_wind_cc_peak.frd`) | **12.78** | n/a (membrane only with `OUTPUT=2D`) | n/a | 0.40 (membrane VM) | | In-house CST + DKT (`tools/shell_dome.py baseline`) | **29.47** | 0.645 | 0.809 | 2.260 | | OpenSees ShellDKGT (`tools/solve_shell_opensees.py baseline`) | **52.42** | 1.745 | 1.179 | 3.007 | Stress recovery for both DKT solvers uses the SAME function (`recover_per_tri_stress`, `src/zomestruct/fea/shell_solver.py:150` — the OpenSees wrapper imports it at `src/zomestruct/fea/opensees_shell.py:63`), so per-fiber σ_b differences mirror displacement differences plus extra amplification from the bending-curvature derivative. Sources: - `reports/full-analysis-shell-opensees-baseline.txt` - `reports/2026-05-04--shell-fea-from-scratch.md` - `reports/2026-05-04--calculix-shell-fea.md` - Re-confirmed by `tools/investigate_shell_side_by_side.py` → `reports/investigate_shell_side_by_side_baseline.txt`. --- ## 2. Step A — flat plate cross-check (already on disk) Both solvers pass the Timoshenko clamped-square-plate reference within the 5% tolerance: | Solver | w_centre / w_theory | Source | |---|---|---| | In-house CST+DKT (`tools/shell_validate_square.py`) | 1.0–1.02 across mesh sizes | `reports/2026-05-04--shell-fea-from-scratch.md` | | OpenSees ShellDKGT (`tools/opensees_shell_smoke_clamped_plate.py`) | 1.024 (10×10 mesh) | `reports/opensees_shell_smoke_clamped_plate.txt` | Element formulation alone is therefore not the bug — both DKT implementations recover the analytical clamped-plate solution. The disagreement is specific to the **curved** geometry. --- ## 3. Step B — per-node disagreement map Tool: `tools/investigate_shell_side_by_side.py baseline`. Output: `reports/investigate_shell_side_by_side_baseline.txt`. Solved both solvers on the SAME mesh + load + BC, computed per-node fractional disagreement `|u_inhouse − u_opensees| / max(|u_inhouse|, ε)`. Classified each node as: - `boundary` — flagged by `MidSurfaceMesh.is_boundary_node` (clamped). - `joint` — multi-panel node (touched by triangles from ≥ 2 panels). - `interior` — intra-panel node. | Class | Count | median |Δu| mm | p95 |Δu| mm | median frac | p95 frac | max frac | |---|---|---|---|---|---|---| | boundary | 338 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | | joint | 911 | 3.22 | 9.22 | 0.264 | 0.555 | 1.027 | | interior | 3 390 | 3.38 | 10.42 | 0.285 | 0.775 | 8.426 | Top-10 nodes by absolute |Δu|: ALL classified `interior`, located in the upper-cap region of the dome at roughly (x, y, z) ≈ (−90 to −600, +850 to +2 080, −1 800 to −2 250). The OpenSees deflection is roughly 3× the in-house deflection at these nodes (e.g. nid 4314: in-house 12.1 mm, OpenSees 41.5 mm). **Diagnosis:** - `boundary` median frac = 0 → BC application is identical (penalty in in-house vs `ops.fix` in OpenSees is not the issue, as expected). - `joint` and `interior` medians are nearly equal (0.264 vs 0.285) → the disagreement is NOT concentrated at panel-joint edges, so the drilling-DOF / inter-element rotation transfer is not the dominant effect. - The pattern is **diffuse** — every interior node disagrees by a similar fraction → the difference is in the per-element formulation itself. This rules out hypothesis (3) — boundary conditions — and weakens hypothesis (1) — drilling-DOF treatment. Step C below confirms (1) is ruled out as the dominant effect. --- ## 4. Step C — drilling-factor sweep Tool: `tools/investigate_drilling_sweep.py`. Output: `reports/investigate_drilling_sweep.txt`. The in-house solver's drilling penalty stiffness is `k_d = drilling_factor * E * t * area` per element, added on the local θ_z DOF only (`src/zomestruct/fea/shell_element.py:367`). The default is `drilling_factor = 1e-2` (`shell_element.py:352`) — note this is NOT the `1e-7` cited in the original task brief; the actual default is 100 000× larger. Sweep results (baseline wind_cc_peak): | drilling_factor | u_max [mm] | σ_b_max [MPa] | VM_max [MPa] | u/u_OPS | u/u_CCX | |---|---|---|---|---|---| | 0 | 30.33 | 23.99 (numerical blow-up) | 27.10 | 0.58 | 2.37 | | 1e-9 | 29.91 | 0.776 | 2.41 | 0.57 | 2.34 | | 1e-7 | 29.91 | 0.776 | 2.41 | 0.57 | 2.34 | | 1e-5 | 29.90 | 0.776 | 2.41 | 0.57 | 2.34 | | 1e-3 | 29.79 | 0.756 | 2.39 | 0.57 | 2.33 | | **1e-2** (current) | **29.47** | **0.645** | **2.26** | **0.56** | **2.30** | | 1e-1 | 28.32 | 0.659 | 1.82 | 0.54 | 2.21 | | 1.0 | 23.88 | 0.629 | 1.40 | 0.46 | 1.87 | Findings: - Across 9 orders of magnitude in drilling stiffness, in-house `u_max` stays within the 24–30 mm band. It NEVER approaches OpenSees's 52.4 mm. - At the extreme `drilling_factor = 1.0` (100× current), in-house only gets STIFFER — moves AWAY from OpenSees, not toward it. - `drilling_factor = 0` triggers numerical instability — σ_b explodes to 24 MPa (singular θ_z mode). The current 1e-2 default is well inside the stable plateau. **Conclusion:** drilling stabilization is not the cause of the in-house ↔ OpenSees disagreement. Hypothesis (1) is ruled out. --- ## 5. Step D — CCX S3 tiebreaker Tool: `tools/investigate_ccx_tiebreaker.py` (parses `reports/ccx_shell/baseline_wind_cc_peak.frd` already on disk; does NOT re-run CCX). Output: `reports/investigate_ccx_tiebreaker.txt`. CalculiX S3 on the SAME mesh, SAME BC, SAME loading: - u_max = **12.78 mm**, u_p99 = 11.58 mm, u_median = 4.56 mm. - Membrane VM_max = 0.40 MPa (CCX `OUTPUT=2D` reports membrane only; through-thickness fiber stress not directly extracted). So all three solvers DISAGREE with each other on `u_max`: | Solver | u_max | factor vs CCX | |---|---|---| | CCX S3 (3-node triangular shell) | 12.78 mm | 1.00× | | In-house CST + DKT | 29.47 mm | 2.31× | | OpenSees ShellDKGT | 52.42 mm | 4.10× | Both pure-DKT solvers over-predict deflection vs CCX by 2–4×. The in-house ↔ OpenSees ratio of 1.78× is part of a larger ~4× spread across the three implementations. --- ## 6. What we suspect — root-cause hypothesis **Confidence:** high on the qualitative diagnosis; medium on the exact quantitative split between contributing factors. The **DKT family of elements** (`Discrete Kirchhoff Triangle`) is a Kirchhoff thin-plate-bending element. Two known limitations matter for faceted curved shells: 1. **No transverse shear stiffness.** DKT is a pure Kirchhoff element — it imposes the Kirchhoff hypothesis (normal-to-mid-surface stays normal) discretely at the corner nodes. On a flat plate, this is the correct simplification. On a faceted approximation of a curved shell, the "missing" shear stiffness across element folds means the discretized shell is artificially soft. 2. **No membrane–bending coupling at the element level.** A DKT element superposes (a) CST membrane + (b) DKT bending + (c) drilling penalty. The three blocks are independent in the local frame. On a curved shell, in-plane stretching IS coupled to out-of-plane curvature change — this coupling is what makes a real shell stiff. DKT only sees this coupling indirectly via the global rotation transformations between non-coplanar adjacent triangles. Different local-frame conventions → different coupling magnitudes. CCX S3, by contrast, is a continuum-based shell with through-thickness shear (per CalculiX's manual it is implemented as a degenerate-S6 with reduced integration — a Mindlin-Reissner rather than Kirchhoff element). On the curved dome it captures the membrane–bending–shear coupling that pure DKT misses, so it predicts a much stiffer response. Why do in-house and OpenSees DKT differ by 1.8× given they're "the same element"? Two suspects: - **Local-frame construction.** In-house uses `e1 = unit(p1 − p0); e3 = unit(cross); e2 = e3 × e1` (`shell_element.py:50-60`). OpenSees `ShellDKGT` uses its own internal local frame; the source is in `SRC/element/shell/ShellDKGT.cpp` in the OpenSees C++ tree, which we don't have at hand. Different local-frame choices → different transformation matrices on non-coplanar joints → different effective stiffness contributions on every shared edge of the curved mesh. Since 19.6% of the dome's nodes are joint nodes (911 of 4 639), small per-edge differences accumulate. - **Drilling-stiffness magnitude.** In-house uses an explicit penalty (1% of membrane scale on θ_z, see Step C). OpenSees `ShellDKGT` source-code reportedly uses zero membrane drilling stiffness (it's the "flat" shell, where rotation about the normal is not assigned a resisting force) — relying on rotational coupling through adjacent non-coplanar elements only. On a flat plate this is irrelevant; on a curved shell it lets the dome rotate more freely at joints and therefore deflect more. This is consistent with OpenSees being SOFTER (more deflection) than in-house, which has a non-zero θ_z penalty. Caveat: the drilling sweep (Step C) shows that even at `drilling_factor = 1e-9` (effectively zero), in-house stays at 30 mm — so removing in-house drilling can't bridge the gap to OpenSees's 52 mm by itself. The drilling difference is therefore at most a contributing factor, not the dominant effect. The dominant effect is most likely **local-frame construction differences interacting with the missing membrane–bending coupling**: on a faceted curved shell, the CST membrane stiffness is computed in a 2D projection that depends on the chosen e1 / e2 axes; rotating that back to global produces stiffness contributions that fight (or assist) the bending block depending on alignment. CCX S3 sidesteps this by using a single curved-shell formulation rather than membrane + bending superposition. --- ## 7. What we ruled out - **(3) BC application.** Step B: median frac disagreement on boundary nodes is 0.000. Penalty (in-house) vs elimination (`ops.fix`) makes no difference here. - **(4) Pressure-load lumping.** Code inspection confirms both compute `area * p / 3 * outward_normal` per-corner translation DOF (`shell_element.py:411-435` for in-house; `opensees_shell.py:139-161` for OpenSees). Identical formula. - **(5) Stress recovery.** Both call the same `recover_per_tri_stress` (`opensees_shell.py:63` imports from `shell_solver.py:150`). Given identical displacement field they produce identical stress. - **(1) Drilling-DOF magnitude (in-house side).** Step C: 9-decade sweep moves in-house `u_max` by < 22%, never approaching OpenSees. Drilling cannot be the dominant cause. The dominant cause is therefore in the **per-element formulation + local-frame construction on the curved mesh** (hypothesis 2 from the brief). --- ## 8. Recommendation for resolution Concrete next steps, in priority order: 1. **Adopt CCX S3 as the primary shell solver for Phase 3 (buckling) and Phase 4 (code-check).** It is the only one of the three with a shell formulation that captures curvature correctly on faceted meshes. The in-house DKT and OpenSees ShellDKGT runs become conservative upper-bound sanity checks. 2. **Switch the OpenSees wrapper to `ShellMITC4` once a centroid-split quad mesher is available** (or use one of the higher-order shell formulations). MITC4 is the OpenSees recommended thin-shell element for curved surfaces — it includes the membrane–bending coupling that DKGT is missing. The current triangular mesh would need a centroid-split into quads (each triangle → 3 quads about its centroid), which is a known follow-up that the existing `opensees_shell.py` docstring already calls out (lines 30-39). 3. **Investigate why OpenSees DKGT is even softer than in-house DKT (1.78× ratio).** Read `SRC/element/shell/ShellDKGT.cpp` from the OpenSees source tree to identify the local-frame convention and confirm whether ShellDKGT applies any drilling stiffness. This is bookkeeping; it does not change the engineering recommendation but would close the loop on the original 2× question. 4. **Re-run the CCX S3 cross-check with `OUTPUT=3D`** to extract peak- fiber stresses (top/bottom surface) so we can compare apples to apples on σ_b, not just `u_max`. Currently CCX reports membrane VM only; estimating fiber stress from membrane requires a 1.5-2× multiplier (see `2026-05-04--calculix-shell-fea.md` § 3). 5. **Mesh-refinement study (Step E, deferred).** Re-run all three solvers at target sizes 100/150/200/300 mm. If CCX S3 converges and the DKT solvers don't, that's the final confirmation that DKT is simply the wrong element for this geometry. Cheap to run; should be done before Phase 3 if budget allows. --- ## 9. Implications for Phase 3 (buckling) and Phase 4 (code-check) The earlier plan (`Option C` from prior planning notes) was to "proceed with both shell solvers' results in parallel" so that Phase 4 has two independent witnesses. Given the findings here, that plan needs revision: - **Buckling (Phase 3).** Buckling factors come from a generalized eigenvalue problem `(K + λ K_g) v = 0` where K is the linear stiffness and K_g is the geometric stiffness from a reference membrane stress state. Since both DKT solvers under-stiffen the shell on the curved dome by 2-4×, their predicted buckling loads will be **non-conservatively under-estimated** (lower buckling load factor than the real shell). Conversely, CCX S3 will give buckling factors closer to physical reality. **Recommendation:** use CCX as the Phase 3 shell-buckling solver. The OpenSees / in-house results may still serve as a "lower-bound" sanity check, but the engineering answer must come from CCX. - **Code-check (Phase 4).** The Phase 4 code-check tabulates per- panel `(σ_b, σ_m, σ_VM)` for ASCE limit-state checks. The two DKT solvers OVER-predict these stresses (by ~2× in-house, ~3-4× OpenSees vs CCX). Using DKT results would be **conservatively unsafe** for the operator (rejecting designs that real shells can carry). For acceptance and final certification, CCX values should drive the D/C ratios. The DKT runs become an "upper bound" that, if it passes, certifies even a more pessimistic interpretation of the shell. - **For "Option C" parallel-track work:** keep the parallel runs only for the validation purpose of demonstrating CCX/DKT spread, NOT for driving Phase 3/4 numerical answers. Document the spread in the acceptance report. --- ## 10. Files Investigation tools created (under `tools/investigate_*.py`): - `tools/investigate_shell_side_by_side.py` — Step B per-node disagreement. - `tools/investigate_drilling_sweep.py` — Step C drilling-factor sweep. - `tools/investigate_ccx_tiebreaker.py` — Step D CCX vs DKT comparison. Captured outputs (under `reports/`): - `reports/investigate_shell_side_by_side_baseline.txt` - `reports/investigate_drilling_sweep.txt` - `reports/investigate_ccx_tiebreaker.txt` External artifacts referenced: - `reports/full-analysis-shell-opensees-baseline.txt` (existing) - `reports/2026-05-04--shell-fea-from-scratch.md` (existing) - `reports/2026-05-04--calculix-shell-fea.md` (existing) - `reports/opensees_shell_smoke_clamped_plate.txt` (existing) - `reports/ccx_shell/baseline_wind_cc_peak.frd` (existing)