Stability issue with P2Form_laplace in 3D
Hi,
I just came across a numeric stability issue with the 3D variant of P2Form_laplace. Apparently for some tetrahedra the code can run into a division-by-zero exception. Below is the output of the HyTeGVsFenicsFormTest when using the following vertex coordinates
vertex | x-coordinate | y-coordinate | z-coordinate |
---|---|---|---|
0 | 0.180901699437495 | 0.131432778029783 | 0.861803398874989 |
1 | 0.180901699437495 | -0.131432778029783 | 0.861803398874989 |
2 | 0.180901699437495 | 0.131432778029783 | 1.111803398874990 |
3 | 0.000000000000000 | 0.000000000000000 | 1.250000000000000 |
and activating floating-point exceptions
[0][INFO ]------(0.018 sec) -----------------------
[0] P2 Laplace Forms (3D)
[0] -----------------------
[0][INFO ]------(0.018 sec) -> running with IDENTITY_MAP
Thread 1 "HytegVsFenicsFo" received signal SIGFPE, Arithmetic exception.
0x000000000051fc49 in hyteg::P2Form_laplace::integrateAll (this=0x7fffffffc930, coords=..., elMat=...)
at /home/mohr/Research/Geophysik/terraneo/Redesign/HyTeG/src/hyteg/forms/form_hyteg_manual/P2FormLaplace.hpp:331
331 INTEGRATE3D( 0, 0 );
The other forms such as P2Form_mass or P2ToP1Form_div do not seem to be affected. Need to figure out how to circumvent this issue.
Cheers
Marcus