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
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
domain), 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
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).
-Jonatan