On 04/29/2015 12:31 PM, Jonatan Lehtonen wrote:
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.
All I've been able to find is the algorithm described in the accepted answer of this StackExchange question <http://scicomp.stackexchange.com/questions/2624/finding-which-triangles-poin...>. 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.
There is also a link to an interesting PDF in Edit2... Anyway, the current code should deal also with hexahedrons with non-planar faces (non-convex!) as long as the centroids stay inside.
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.
Yes, I have thought about this possibility, but decided not to go that way now :).
The best option is probably to go with what you've got and add the fall back brute force, as you described.
Or using quad-/oct-trees with the brute force (but kdtree is available and fast).
Try the new code [1], please, and let me know if it works for your (non-awful) meshes. I have also fixed the nan vs. zero problem.
The brute force will have to wait.
r. [1] https://github.com/rc/sfepy/tree/improve-find-ref-coors