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.
>>
>> [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
> <http://scicomp.stackexchange.com/questions/2624/finding-which-triangles-points-are-in>.
> 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