Heat_method_3: fix multi-source distance normalization (value_at_source_set)

[jf---]

Surface_mesh_geodesic_distances_3::estimate_geodesic_distances() returns badly wrong distances for multi-vertex source sets — the more sources (and the larger the spread of the recovered potential across them), the worse. The flagship "distance from the domain boundary" usage (add every boundary vertex of an open mesh as a source) reads ~0.5x the true distance near the boundary, with most near-boundary vertices reading below the euclidean straight-line lower bound, which is impossible for a distance. Single-source results are correct, which is why the examples and the existing tests never exposed it.

Root cause. The Poisson step recovers the potential phi only up to one additive constant, and phi is in general not constant across a multi-vertex source set (the normalized gradient field is not exactly integrable). value_at_source_set() then computes

d(i) = min over sources s of |phi(i) - phi(s)|

which is the 1-D distance from the value phi(i) to the set of source values, not a shift: every vertex whose phi lands inside the source-value spread is folded toward 0, and the far field is shifted by max_s phi(s) instead of a single constant. With one source the rule degenerates to the correct shift.

Fix. Anchor phi with a single constant: the mean of phi over the source set (the least-squares choice, matching common practice in other heat-method implementations, e.g. geometry-central). Single-source behavior is unchanged (mean of one value = that value). Distances at source vertices are now ~0 up to discretization error instead of exactly 0-by-construction of the fold; the existing package tests (which check distance ranges) pass unchanged.

Test. heat_method_surface_mesh_multisource_test.cpp + data/flat_star_disk_irregular_rim.off: a small flat star-shaped disk (z = 0) with irregular, variable-density boundary sampling, all 508 boundary vertices used as sources. Because the mesh is flat and simply covered, the exact geodesic distance to the source set equals — in the near-boundary band — the euclidean distance to the nearest source (verified to ~1e-15 against exact polyhedral geodesics when the fixture was generated; agreement stated in the file header), and the euclidean distance is a hard lower bound everywhere. So the test is self-contained: no oracle dependency, and the failure is unambiguous (a distance smaller than the straight line). The fixture is fully synthetic (generated programmatically; CC0).

Measured with this test (band = 2h..6h of the euclidean distance to the source set, h = mean edge length), identically in Direct and Intrinsic_Delaunay modes:

	band mean d / d_true	band vertices below the euclidean lower bound
before (min over source values)	0.544	83.4 % (assertions fail)
after (mean shift)	1.135	2.4 % (marginal)

On a real-world case (12.9k-vertex planar slice of a Stanford-bunny print, 434 boundary sources, measured against exact MMP geodesics): near-boundary ratio 0.41 before, 1.00 after.

[sloriot]

Successfully tested in CGAL-6.3-Ic-31