[BUG]: GJK fails on simple convex hull

[jorgensd]

Summarize the issue

Running the following minimal example makes GJK not converge.
First reported at:
https://fenicsproject.discourse.group/t/gjk-error-in-interpolation-between-non-matching-second-ordered-3d-meshes/16086/8?u=dokken

How to reproduce the bug

See below

Minimal Example (Python)

import numpy as np


p0 = np.array([[44.2919, 39.7508, 18.8208]])
p1 = np.array([[46.0659, 38.6251, 16.3637], 
                 [44.2764, 42.3867, 18.239],
                 [  50.0397, 38.4649, 19.9843], 
                  [ 47.1963, 40.4563, 19.1057] ,
                  [ 48.0266, 38.5593, 18.2475] ,
                   [45.1735, 40.495, 17.3336] ,
                   [48.0415, 40.4456, 18.1963 ]])
print(p1.shape)
import dolfinx
from dolfinx.geometry import compute_distance_gjk

compute_distance_gjk(p0, p1)

Output (Python)

Traceback (most recent call last):
  File "/home/dokken/Documents/src/debug/mwe.py", line 16, in <module>
    compute_distance_gjk(p0, p1)
    ~~~~~~~~~~~~~~~~~~~~^^^^^^^^
  File "/home/dokken/src/mambaforge/envs/test_mpl/lib/python3.13/site-packages/dolfinx/geometry.py", line 272, in compute_distance_gjk
    return _cpp.geometry.compute_distance_gjk(p, q)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^
RuntimeError: GJK error - max iteration limit reached

Version

main branch

DOLFINx git commit

No response

Installation

Tested both docker and conda, stable and nightly

Additional information

No response

[jorgensd]

Further reduced to:

import numpy as np


# p0 = np.array([[44.2919, 39.7508, 18.8208]])
p0 = np.array([[44.2919, 39.7508, 18.8208]])
p1 = np.array(
    [
        [46.0659, 38.6251, 16.3637],
        [44.2764, 42.3867, 18.239],
        [50.0397, 38.4649, 19.9843],
        [45.1735, 40.495, 17.3336],
        [48.0266, 38.5593, 18.2475],
    ]
)
import matplotlib.pyplot as plt
import matplotlib

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection="3d")
ax.scatter(p0[:, 0], p0[:, 1], p0[:, 2], c="r", marker="s", label="p0")
ax.scatter(p1[:, 0], p1[:, 1], p1[:, 2], c="b", marker="o", label="p1")
ax.legend()
for i in range(p1.shape[0]):
    ax.text(p1[i, 0], p1[i, 1], p1[i, 2], f"p{i}")
plt.savefig("test.png")


print(p1.shape)
import dolfinx
from dolfinx.geometry import compute_distance_gjk

compute_distance_gjk(p0, p1)