geoid harmonic projection re-derives the surface mesh with ConvexHull instead of integrating over the boundary facets we already have

[lmoresi]

geoid._project_spherical_harmonic_samples (postprocessing/geoid.py:240) computes the
harmonic coefficient by gathering every boundary coordinate and value to rank zero,
discarding the mesh connectivity, and re-deriving a triangulation of the point cloud
with scipy.spatial.ConvexHull(unit_coords, qhull_options="QJ"). It then integrates with
analytic spherical triangle areas and one-point centroid values.

The boundary facets already exist — it is the FE trace, and _desmear walks exactly
those facets a few modules away. uw.maths.BdIntegral already integrates over them
properly, through DMPlexComputeBdIntegral, with the real quadrature and no gather.

The coefficient is then just the projection onto the harmonic written in the FE inner
product the assembler already uses:

numerator   = uw.maths.BdIntegral(mesh=mesh, fn=h.sym[0] * legendre, boundary=b).evaluate()
denominator = uw.maths.BdIntegral(mesh=mesh, fn=legendre**2,        boundary=b).evaluate()
coefficient = sign * numerator / denominator

What is wrong with the current path

• It cannot work on a non-convex boundary. A convex hull of the boundary nodes is
  only the boundary for a convex surface. An annulus arc, an internal boundary, or a
  deformed free surface will silently get the wrong triangulation.
• QJ joggles the input — it perturbs coordinates to break degeneracies. That is a
  poor foundation for a reproducible diagnostic.
• It mixes measures. Analytic spherical triangle areas are applied to a faceted FE
  reaction; the two disagree at O(h²), which is the chord-versus-arc family of Curved-boundary P2 VERTEX flux/traction values converge only slowly — document and steer consumers to midpoint/fitted quantities #414/CBF sigma_nn recovery on deformed/faceted boundaries: ~2.4% amplitude deficit (chord-vs-arc boundary mass suspected) #431.
• It gathers the whole boundary to rank zero and triangulates there — the same
  scaling objection as boundary_flux 3D assembly gathers the full global boundary onto every rank (memory scaling) #408.
• Centroid value = mean of three nodal values is a one-point quadrature.

Measured

Zhong l=2 shell, same solves, coefficient three ways, relative error against the analytic
0.41920 (surface) / 0.77060 (CMB):

cellSize	0.25	0.20	0.16	0.13
surface, ConvexHull	0.0028	0.0001	0.0014	0.0013
surface, BdIntegral	0.0003	0.0014	0.0005	0.0008
CMB, ConvexHull	0.0041	0.0081	0.0059	0.0027
CMB, BdIntegral	0.0091	0.0098	0.0084	0.0037
CMB, pointwise fit	0.0313	0.0236	0.0188	0.0097

On accuracy the two integral arms are a wash — ConvexHull is slightly better at the
CMB, BdIntegral slightly better at the surface, neither converging cleanly. This is not
an accuracy argument and should not be sold as one. It is a correctness-of-construction
and portability argument: same answer, on the mesh we actually have, on any boundary
shape, without a gather or a joggle.

The third row is worth noting separately: both integral forms beat the pointwise harmonic
fit by 3-10x at coarse resolution, which is independent support for #414's advice to
steer consumers to integral/fitted quantities rather than sampled vertices.

Relation to #646

#646 adds boundary_normal_traction_integral, which contracts the assembled reaction
with the test function directly and needs no pointwise recovery at all. For a coefficient
that is better than either arm above, because it never forms a pointwise field. This
issue is about the path that remains when projection="centroid" — currently the
default. Between the two, _project_spherical_harmonic_samples looks retirable rather
than repairable.

Script: ~/+Simulations/topography_penalty_633/replace_convexhull.py

[lmoresi]

Preferred resolution (Louis, 2026-08-26): use UW machinery, and @gthyagi's reaction projection in #646 is the one to take. Recording that here so this issue does not read as advocating my BdIntegral sketch over it.

The ordering is:

1. @gthyagi's boundary_normal_traction_integral (Add distributed harmonic projection of rotated reactions #646) — contracts the assembled reaction with the test function directly, over owned DOFs plus MPI reductions. For a harmonic coefficient this is strictly the best of the three: it never forms a pointwise field, so the zero-mean P2 vertex basis (Curved-boundary P2 VERTEX flux/traction values converge only slowly — document and steer consumers to midpoint/fitted quantities #414), the singular vertex rows of the P2 triangle mass (Grad-div penalty breaks POINTWISE traction recovery (28% on dynamic topography); facet-integrated is fine #633) and the chord-versus-arc geometry error (CBF sigma_nn recovery on deformed/faceted boundaries: ~2.4% amplitude deficit (chord-vs-arc boundary mass suspected) #431) are all simply not in its path. It is also already distributed and tested.
2. uw.maths.BdIntegral — the UW-native boundary integral, for the cases where a pointwise field genuinely exists and has to be integrated (the free-surface driver needs a nodal field regardless). Same argument as above about using the mesh we have rather than re-deriving one.
3. _project_spherical_harmonic_samples — retire. Nothing it does is unavailable from (1) or (2), and it carries the gather, the QJ joggle, the convexity assumption and the analytic-sphere-versus-faceted-FE measure mismatch.

The measurements in the issue body stand, with the caveat already stated there: between (2) and (3) the accuracy is a wash, so the argument was never accuracy. It is that (1) and (2) are UW machinery — parallel, tested, and exercised across the codebase — while (3) is a bespoke rank-zero geometry routine that only works on a convex surface.

Suggested action: fold the retirement into #646 if that is convenient for @gthyagi, otherwise a follow-up once #646 lands.

Underworld development team with AI support from Claude Code