fix: correct vertex singularity detection in computeSingularityTerms

[sasso-effe]

Summary

• Floating-point arithmetic could cause computeSingularityTerms to select the edge singularity path (Case 2) instead of the vertex singularity path (Case 3) for query points whose projection P' lands at a face vertex.
• The g1/g2 branch selection in Case 3 used an exact == 0.0 comparison inconsistent with the < EPSILON_ZERO_OFFSET threshold used everywhere else in the function.

Root cause

When P' projects onto a vertex, floating-point arithmetic produces rvn[(j+1)%3] ≈ ~|G| - eps instead of exactly |G|. This satisfies the Case 2 condition rvn < |G|, firing the edge singularity formula -pi * h_p instead of the correct vertex formula -theta * h_p.

Changes

• Case 2 guard: added >= EPSILON_ZERO_OFFSET on both vertex-distance norms so near-vertex points (norm < ε) fall through to Case 3.
• Case 3 g1/g2 selection: replaced r1Norm == 0.0 with r1Norm < EPSILON_ZERO_OFFSET for consistency with the rest of the function.

Verification

Create a right-angle tetrahedron with vertices (0,0,0), (2,0,0), (0,2,0), (0,0,2) and query point (10, 0, z).
Changing z and comparing the potential produced by original vs fixed you obtain:

z	original	fixed
1e-10	9.34e-12	9.34e-12
1e-12	9.34e-12	9.34e-12
1e-13	9.34e-12	9.34e-12
1e-14	-5.23e-09	-5.23e-09
1e-15	-5.23e-09	9.34e-12
sin(pi) ≈ 1.22e-16	-5.23e-09	9.34e-12
1e-16	9.34e-12	9.34e-12
1e-17	9.34e-12	9.34e-12
0.0 (exact vertex)	9.34e-12	9.34e-12

Discontinuities are highlighted in bold.

All existing tests pass.

Remaining bugs!!!

z = 1e-14 shows a pre-existing discontinuity in both the original and fixed code, which I was not able to fix. It could be related with EPSILON_ZERO_OFFSET being exactly 1e-14

[schuhmaj]

Thanks for reaching out and opening a PR.

The reason that 1e-14 is still off lies in the fact that the sum2 containing the transcendental $AN_{pq}$ is not inline with the evaluation of $sing_{\mathcal{A}}$.

The extra condition for case 2 is definitely a good addition. For the floating point comparison, I am thinking of inroducing some sort of almostZero function. Though introducing such a function led in a first test run to consistent results, i.e., fixing the discontinuity, but different values - i.e., the tests failed.

Otherwise, the question is, anyway, how much we can fix this. Because with any EPSILON, we will always have some sort of discontinuity somewhere.
So far, the strategy was to throw a warning if the floating point magnitudes are becoming too large, and elimination can become a thing (in the isCriticalDifference function which happens from the 1e-16 case onwards since, we use a difference in exponents of 50, i.e. $2^{-50} ≈ 8 \times 10^{-16}$, which is not in-line with our other 1e-14 epsilon at the moment). So another fix could also be to just recalibrate when this warning triggers.
(Ofc, the the changes you made here would then still be good to include 😀)

Anyways, I would leave this PR open for now until we have a full fix and will continue here next week (have other todos around my ears as well unfortunatley)

---

Tetrahedron Potential Evaluation — Before vs. Fix

Each block lists the calculation point P and the resulting potential, for the
original run (before) and the run after the absolute-epsilon fix (fix).
For each of the tetrahedron's four planes, the values sum2, planeDistance and sing
are reported. Columns are merged (colspan) where both runs agree; split where they differ.

P = [10, 0, 1e-13]

potential = 9.344421547810072e-12 (before = fix)

Plane	sum2		planeDistance		sing
	before	fix	before	fix	before	fix
1	-2.47549e-16		1e-13		0
2	0		0		0
3	3.12198		10		-31.4159
4	-0.0196091		4.6188		0

P = [10, 0, 1e-14]

potential = -5.2326385404157645e-09 (before = fix)

Plane	sum2		planeDistance		sing
	before	fix	before	fix	before	fix
1	-2.47549e-17		1e-14		0
2	0		0		0
3	⚠️1.55119		10		⚠️-31.4159
4	-0.0196091		4.6188		0

P = [10, 0, 1e-15]

potential: before = -5.2326385404157786e-09, fix = 9.344421547809656e-12

Plane	sum2		planeDistance		sing
	before	fix	before	fix	before	fix
1	0		1e-15		0
2	0		0		0
3	1.55119		10		-31.4159	✅-15.708
4	-0.0196091		4.6188		0

P = [10, 0, 1e-16]

potential = 9.344421547807886e-12 (before = fix)

Plane	sum2		planeDistance		sing
	before	fix	before	fix	before	fix
1	0		1e-16		0
2	0		0		0
3	1.55119		10		-15.708
4	-0.0196091		4.6188		0

P = [10, 0, 0]

potential = 9.344421547807886e-12 (before = fix)

Plane	sum2		planeDistance		sing
	before	fix	before	fix	before	fix
1	0		0		0
2	0		0		0
3	1.55119		10		-15.708
4	-0.0196091		4.6188		0

[darioizzo]

@schuhmaj Merging this as we need the "temp" fix to have a cleaner workflow for #54 . Its a small and ineffective change anyway, you can still try to provide a more generic solution when you have time.