Yes, it's horrible. :)
Short story: The new code fails as well, I know why (non-convex domain), and do
not know a workaround (quick internet search did not help either).
Long story: Imagine a mesh with a hole (e.g. a long rectangular hole). The mesh
has nice small elements on one side of the hole, and very big elements on the
other side, much bigger than the hole side-to-side dimension. If a target
(probe) point is in a big element, and the closest mesh point in a small
element on the other side of the hole, I have no idea, how to reach it from
there. The current algorithm follows the elements in the direction from the
closest point to the target, but it cannot cross the hole. It would work for a
convex hull mesh of the original mesh, though.
So, for the moment, I will merge the current "solution", as it is about 2-3x
faster than the old one and works well for reasonable meshes (much better than
the old one).
We may also add a slow "brute force" fall back algorithm (search all elements
in spheres of increasing radius centered in the target point - scipy's kdtree
will help here), that will be optionally tried for points not found by the fast
one.
If you are interested, try searching for articles or for how libraries like VTK
do it (there are probes in paraview, and VTK probes are now in sfepy as well,
thanks to Vladimir, see [1], but work only for approximation order 1 (or 0)).
r.
[1] http://sfepy.org/doc-devel/primer.html#probing
All I've been able to find is the algorithm described in the accepted answer of
this StackExchange question. The upside is that it should never fail for awful.mesh. The downside is that it may fail to converge in non-convex meshes (since it cannot jump over gaps), it might only work for triangular/tetrahedral meshes, and it runs in O(N^(1/d)). The runtime would probably be significantly improved by first finding a nearby mesh node with the kdtree, though.
My intuition says that it's always possible to find a counter-example for any non-brute-force algorithm though, so it may not be worth spending too much time on it. E.g. while the above algorithm would work for awful.mesh, it would fail the moment you add one more triangle with corners at (0, 1), (-0.1, 4) and, say, (0.5, 0.6). Then it wouldn't converge for a point at (0.45, 0.45) with first-order elements.
Naturally one option would be to add an extra hidden mesh which, when combined with the original mesh, gets the convex hull of the object. That way you could probably guarantee convergence, but that would involve a considerable computational cost.
The best option is probably to go with what you've got and add the fall back brute force, as you described.
-Jonatan