Points on a convex curved boundary are always extrapolated, silently, with O(h) error

[lmoresi]

Every query point on a convex curved boundary is extrapolated rather than
located, at every rank count including serial, and the caller is not told.

This came out of the #604 / #602 discussion about whether points on the domain
boundary are findable. The answer is that they are not, and that it is not a
parallel defect.

Measured

Annulus, radiusInner=0.5, radiusOuter=1.0, cellSize=0.15. A P2
MeshVariable carrying x^2 + 2 y^2, which P2 represents exactly, sampled at 36
angles chosen to fall BETWEEN boundary nodes. check_extrapolated=True reports
per point whether the value was located or extrapolated.

ring	extrapolated	max abs error
outer, r = R exactly	36/36	5.97e-02
outer, r = R(1 - 1e-9)	36/36	5.97e-02
outer, r = R(1 + 1e-9)	36/36	5.97e-02
inner, r = R exactly	0/36	2.8e-16
mid-annulus	0/36	6.7e-16

Identical at np=1, np=2 and np=4 — the numbers above do not change with rank
count.

Why, and why the two boundaries differ

The mesh is a polygon. On the OUTER boundary its chords sag inside the true
circle, so a point on the true circle between two boundary nodes is outside the
meshed domain by about the sagitta. On the INNER boundary the chords lie closer
to the centre than the true circle, so the meshed domain extends slightly further
in and a point at r = R_in is inside it. Which side the domain is on decides
whether "on the boundary" means inside or outside.

A point a nanometre inside the nominal outer radius is still outside the mesh,
which is why tightening the query radius does not help.

The error converges at O(h), not O(sagitta)

cellSize	extrapolated	max abs error	sagitta ~ h^2/8
0.300	36/36	1.095e-01	1.13e-02
0.150	36/36	5.970e-02	2.81e-03
0.075	36/36	2.668e-02	7.03e-04

The miss distance falls as h^2 but the error falls as h, so this is the
extrapolant's own error over a stencil of size h rather than the distance
travelled. Refinement helps, slowly: still 2.7% at cellSize=0.075.

Why it matters

evaluate and global_evaluate return return_value alone unless the caller
asks for check_extrapolated=True. So sampling a field on the outer boundary of
an annulus or a spherical shell — surface velocity, dynamic topography, a
boundary profile for a benchmark — silently returns an extrapolated value with a
first-order error, and nothing in the call says so.

Directly relevant to #320 (free-slip methods compared on Kramer annulus and
spherical benchmarks) and to anything reading a field on the outer surface of a
shell.

What the fix might be

Not obvious, and it is a contract decision rather than a bug fix:

1. Snap query points to the mesh. Pull a point that misses by less than a
   cell inward to the nearest facet, and interpolate there. Cheap, and it makes
   r = R mean "on the meshed boundary" rather than "on the ideal circle".
2. Say so. Warn once per call with a count when any point was extrapolated.
   Does not fix the number, but it stops the error being invisible — and it is
   the same request as A located point can interpolate to NaN, and #602 makes that invisible #604's.
3. Leave it and document it, with the guidance that boundary sampling should
   use the mesh's own boundary nodes.

Option 2 is worth doing whatever else is decided.

Not a parallel problem

Worth stating plainly, because that is where the suspicion started: the parallel
best-claim fallback is working. Serial and parallel agree to the digit on every
row above, so the machinery that recovers unlocatable points is doing its job.
What is wrong is that a point on a curved boundary is genuinely outside the
domain we meshed.

Related: #604 (a located point interpolating to NaN — a different mechanism),
#551 (the locator umbrella), #320.

Underworld development team with AI support from Claude Code

[lmoresi]

Closing this — it was raised by us and it does not stand up as an issue.

The core finding is not a defect. A linear mesh is a polygon, so a point on the true circle between two boundary nodes is outside the domain that was actually meshed, and extrapolating it is the only available answer. There is nothing here to fix short of giving boundary cells a curved (higher-order) coordinate space, which is a change of discretisation rather than a bug fix.

The measurements stand and are worth keeping visible, so restating them here rather than losing them with the issue: on an annulus at cellSize=0.15, 36 of 36 points at r = R_outer are extrapolated at np=1, 2 and 4, with a max error of 5.97e-02 on a quadratic P2 field; points at r = R_outer(1 - 1e-9) behave identically; the inner radius is located exactly (0/36, 2.8e-16) because there the chords lie inside the true circle and the meshed domain extends past it. The error falls as O(h) — 1.09e-01, 5.97e-02, 2.67e-02 at h = 0.30, 0.15, 0.075 — rather than with the O(h^2) sagitta, so it is the extrapolant's own error over an h-sized stencil and refinement helps only slowly.

The single actionable residue is that the extrapolation is silent: neither evaluate nor global_evaluate returns the mask unless asked. That is already requested in #604 and #606, so it does not need an issue of its own.

The practical guidance — sample on the mesh's own boundary nodes rather than on the ideal surface, and expect an extrapolated value if you sample the ideal one — belongs in the documentation, not in a tracker.