This is quite suprising, as this indicates that the point is not within the cell, as the reference coordinate should be within the reference tetrahedron.
This is due to the fact that you only compute the distance to the bounding box, and not the convex hull of the cell, i.e.
should be
candidates = geometry.compute_collisions_points(tree, pt[None, :])
candidate_cells = geometry.compute_colliding_cells(m_src, candidates, pt.reshape(-1, 3))
candidate_cells = candidate_cells.array[candidate_cells.offsets[0]:candidate_cells.offsets[1]]
Furthermore, by inspecting the distance to the function in reference space, you get:
best_val = 0.0
max_ref_distance = 0.01
import basix
ref_geom = basix.geometry(m_src.basix_cell())
# Iterate through candidates to find the most geometrically valid cell
for c in candidate_cells:
print(geometry.compute_distance_gjk(m_src.geometry.x[m_src.geometry.dofmap[c]], pt))
g_idx = mesh.entities_to_geometry(m_src, tdim, np.array([c], dtype=np.int32))[0]
coords = m_src.geometry.x[g_idx]
try:
# Pull back to reference space to check for extrapolation violation
ref = cmap.pull_back(pt[None, :], coords)[0]
dist = np.linalg.norm(geometry.compute_distance_gjk(ref, ref_geom))
if dist > max_ref_distance:
breakpoint()
as you can observe, the node lies within the convex hull of the cell, but not the cell itself.
It should lie within another cell:
I’ve made a PR to fix this: Improve non-matching interpolation for non-affine grids. by jorgensd · Pull Request #4127 · FEniCS/dolfinx · GitHub
where the undershoot reduces from -0.47 to -0.09, which is more expected as you are interpolating data that might be negative with a cell.