On Wednesday, April 29, 2015 at 3:18:59 PM UTC+3, Robert Cimrman wrote:
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.

[1] https://github.com/rc/sfepy/tree/improve-find-ref-coors

I tried your solution with the mesh I've been working on, using roughly 5 million measurement points inside the mesh. With the old code (the latest major release, I think), probing took 46 seconds and resulted in 898 NaNs. With the new code, it took 20 seconds and resulted in no NaNs. So it certainly looks like your code works. =)

While working with the new code, I noticed that apparently get_evaluate_cache has been removed from FEDomain. Currently in my code, after solving the Laplace problems with code very similar to the OP, I execute the following line of code:

Probe.cache = pb.domain.get_evaluate_cache(cache=None, share_geometry=True)

So what I'm wondering is, what should I replace this code with? Essentially what I want to do is ensure that every probe I create has the evaluate cache corresponding to that mesh (so that if I were to run the code again with a different mesh, every subsequent probe would automatically use the evaluate cache of the new mesh).