# Geometry failure diagnosis — 2026-09-08 ## Conclusion The first production Newton correction is small in surface amplitude but produces large, oscillatory radial displacement gradients near the eight core corners. The tensor-P3 representation of the prescribed extension distorts the intended ray-wise profile; the poorly conditioned reference mesh amplifies its gradients. There is also a verified geometry-safety defect: the production quadrature union misses inversion between its samples. At the first accepted step length `alpha = 0.086127771767489855`, the relative mapping determinant is **-0.670218 at the core corner and -0.428915 at an interior point close to it**, while the reported limiting quadrature point has determinant **+0.100006**. Thus this is not merely an approach toward a future singularity: the first accepted step already corresponds to a folded represented geometry. The production operators can continue evaluating because their own sampled points remain admissible. No physical equations, mapping/extension implementations, Newton policies, or production tolerances were changed. This investigation implemented an opt-in experiment and an internal, guarded diagnostic-access hook only. ## Reproduction and artifacts Build: ```sh cmake --build cmake-build-release-homebrew --target geometry_quality_experiment -j 6 ``` Runs performed, all single-rank Release: ```sh ./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_baseline --tolerances 0.03,0.003 --max-linear-iterations 100 --save-vectors ./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_profile --replay-vectors geometry_quality_results_2026-09-08_baseline/newton_0_vectors.csv --no-fd --no-block-actions ./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_vertices --replay-vectors geometry_quality_results_2026-09-08_baseline/newton_0_vectors.csv --diagonal-only python3 experiments/summarize_geometry_quality.py geometry_quality_results_2026-09-08_baseline python3 experiments/summarize_geometry_quality.py geometry_quality_results_2026-09-08_vertices ``` The runs took 510.721, 152.037, and 92.3526 seconds, respectively. Replay checks the saved accepted state against a fresh context and verifies `Jp+F`, rather than re-solving for a new correction. Output directories are refused if they already exist. Help, invalid-tolerance rejection, existing-directory rejection, and normal context restoration were exercised. `git diff --check` passed. No test suite was built or run. The existing build cache contained a nonexistent Eigen 5.0.1 include path. CMake dependency discovery was refreshed with `-U '*eigen3*'`, locating installed Eigen 3.4.0. See the run provenance and metadata files. The fresh baseline reproduces the earlier sandbox's seed residual, 23 Krylov iterations, and first safe step. ## 1. Element 0 is the representative of a core-corner pattern The eight limiting elements are **0, 9, 18, 27, 36, 45, 54, 63**. They are core elements (attribute 1), not compactified exterior elements. For the Newton direction their quadrature boundaries differ by only about 0.013%; element 0 has the smallest value. They are near-ties, not exact ties. For uniform contraction they agree within one part per million. The actual Newton limiting sample is on element 0's body diagonal: - Integration coordinates: `(0.010885670927, 0.010885670927, 0.010885670927)`. - Physical reference position: approximately `(-0.144337126, -0.144337126, -0.144337126)`. - Physical radius: `0.249999235`, immediately inside the core boundary. - Reference element Jacobian singular values: minimum `5.51738e-5`, maximum `0.108818`, condition number approximately **1,972**. - At the exact corner the corresponding condition number is approximately **3,151**. The relative displacement mapping at the seed is the identity. Its determinant of one therefore conceals the conditioning of the underlying reference element. Evidence: baseline `newton_0_geometry_elements.csv`, `uniform_contraction_reference_corner_probes.csv`. ## 2. The dangerous correction is small and nonuniform Unweighted surface-parameter statistics for the first correction: | Quantity | Fraction of reference radius | | --- | ---: | | Mean | -0.000201759 | | RMS | 0.000360659 | | Nonmean RMS | 0.000298945 | | Most negative | -0.000751271 | | Most positive | +0.000661834 | The largest displacement is only **0.0751% of the radius**. The negative extrema occur at the eight cube-corner surface directions; positive extrema occur near the twelve edge-center directions. The pattern is nearly invariant under cube symmetries: the largest spread among symmetry-equivalent nodal directions is `1.3544e-6`, compared with the full amplitude range `0.0014131`. The mean-only correction has no sampled boundary through alpha=1. Removing the mean changes the quadrature boundary from `0.0956975242` to `0.0957208926`, only 0.0244%. The nonmean component therefore accounts for almost the entire local compression. It is not an oversized uniform radius change or a warm-start-only artifact; the first solve starts cold. ## 3. Direct diagonal inspection identifies the interpolation mechanism On the negative core-corner ray, the target surface point is fixed. For the default radial-power-2 prescription and logical stellar radius one, the intended continuous radial displacement is exactly `u_r(s) = a_corner * r_logical(s)^2`, with `a_corner = -0.0007504485104433`. At the Newton limiting sample: | Quantity | Intended continuous ray profile | Represented P3 FE field | | --- | ---: | ---: | | Radial displacement | -4.6393850e-5 | -3.5188044e-5 | | Physical radial derivative | -0.4881315 | **-10.4495911** | The derivative is amplified **21.4 times relative to the intended profile**. The values agree at the P3 interpolation nodes, but disagree between them. At `s=0.1`, the FE radial displacement even becomes positive although the intended profile remains negative. The code constructs the extension from logical-radius weights and surface-trace interpolation at volume nodes, then represents it in the tensor-P3 volume space. The corner element spans three logical max-coordinate sectors, projecting toward three surface faces. The resulting angular/radial composition is not a single low-degree tensor polynomial inside that element. Its nodal interpolant need not preserve the prescribed ray-wise radial behavior. The physical gradient further multiplies by the inverse reference Jacobian. Near the core corner, the physical radial coordinate changes extremely slowly along the logical ray. At the limiting sample, `dr_logical/ds=-0.125` but `dr_physical/ds=-9.55639e-5`. This explanation is supported by the measured saved mesh, not just by mesher source. Sibling STROID source suggests why the core map flattens there, but its working tree is dirty and was not treated as proof of mesh provenance. Evidence: profile/vertices `newton_0_core_diagonal.csv`; production implementation `libmeanfield/impl/deformation/radial_extensions.cpp` and field registry P3 displacement versus saved mesh P4 geometry. ## 4. Quadrature safety is not element safety The unchanged production estimator was run with a second sampling plan containing the eight vertices of each of the eight core-corner elements. Independently, the diagnostic formed the polynomial `det(J_reference + alpha * dU/dxi) / det(J_reference)` directly from element matrices, without inverse Jacobians. The methods agree. | Direction | Production quadrature boundary | Vertex boundary | | --- | ---: | ---: | | Original Newton correction, radial power 2 | 0.0956975 | **0.0515682** | | Same surface correction, radial power 3 | 0.382625 | **0.206273** | | Same surface correction, radial power 4 | No boundary through 1 | **0.825091** | For the original correction at alpha=0.0861277718: | Diagonal coordinate s | Relative determinant | | --- | ---: | | 0, exact corner | **-0.670218** | | 0.005, interior point | **-0.428915** | | 0.01088567, production limiting quadrature point | +0.100006 | | 0.02 | +0.609367 | The negative interior value establishes an actual fold, not merely a singular boundary vertex. Adding vertices would catch this case, but vertex sampling alone is not a general positivity certificate: the uniform-contraction control has a much stricter interior quadrature boundary than vertex boundary. ## 5. Controls and ruled-down explanations ### Linear accuracy | Requested linear tolerance | Krylov iterations | Verified relative residual | Geometry-safe step | | --- | ---: | ---: | ---: | | 0.03 | 23 | 0.0244757 | 0.0861278 | | 0.003 | 37 | 0.00299622 | 0.0865524 | Tenfold tightening changes the normalized surface correction by 1.30% and the safe step by only 0.49%. Linear inaccuracy is not the main explanation for this first-step failure, although this does not establish behavior at every later state or arbitrary tolerance. ### Mapping and derivative consistency - Fresh mapped geometry agrees with the preflight affine prediction to about `6e-16` for the actual Newton correction. - Finite differences of generated volume displacement agree with the production extension direction to about `1.6e-16`. - Material/shape/hydrostatic directional residual differences improve roughly tenfold with the first tenfold perturbation reduction, reaching relative errors around `1e-7`, then show roundoff/cancellation. - Gravity checks improve initially but have smaller net Jacobian actions and noisier relative errors. Mass finite differences are cancellation-sensitive; the smallest perturbation is worse. These checks do not prove the entire Jacobian correct, but reveal no gross inconsistency along this correction. - The central-density relative error of one has absolute magnitude `1.22e-23`; it is not evidence of an important phase-border defect. ### Extension and representation controls Uniform contraction through the default extension has a quadrature boundary of 0.04819 (4.82% surface contraction), independent of Newton/EOS balances. Core-only P3 projections of physical affine and physical-radius-squared displacements have boundaries 0.2014 and 0.7693. These are diagnostic bypass controls, not production solutions. Even physical affine scaling is not exactly representable by P3 displacement on P4 geometry. Higher radial powers reduce core sensitivity, but **power 4 is not a demonstrated fix**: despite all production quadrature points permitting alpha=1, its corner determinant at alpha=1 is **-0.211997**, and a nearby interior determinant is also negative. Increasing order or power alone should not be assumed to solve the problem. In particular, faithfully interpolating the original logical-radial profile can make the uniform control more singular near a poorly conditioned reference corner; interpolation sometimes suppresses as well as amplifies it. ## Recommended next work 1. Add this replay direction and the negative interior/corner determinant as a targeted geometry regression. Extend admissibility sampling to vertices and near-corner/edge regions, then consider adaptive or bounded determinant certification for high-order elements. Do not simply reduce the minimum step. 2. Design an extension that controls gradients in physical reference coordinates and preserves simple radial/affine modes in the represented FE space. Treat reference-mesh quality and displacement/geometry order compatibility together. 3. Rebuild the coupled residual/Jacobian and solve anew under any changed extension. Do not accept a post-processed old Newton direction as though it solved the new Newton equation. 4. Recheck the mesh-symmetric nonuniform surface pattern and block residuals after repairing the geometry representation. The precise origin of that discrete pattern, and full nonlinear convergence after a repair, remain unestablished. The immediate geometric failure mechanism and the sampling defect are now reproduced. A production remedy and its nonlinear convergence are deliberately not claimed by this diagnostic experiment.