On 04/28/2015 07:52 PM, Jonatan Lehtonen wrote:

Hi Robert,

Sounds great!

Unfortunately I don't really know how to manually create arbitrary meshes in gmsh (the meshes it generates by default are "too good" for this purpose). However, here's at least one mesh that I think should be fairly horrible (and is aptly named). There are some points in there for which there are almost ten elements "between" the point and its closest mesh node. Naturally, you can make the mesh worse by an arbitrary factor by changing nodes 4-13 (increase y-coordinate or push x-coordinates closer to each other) or by adding more triangles in the same style.

I tested this mesh with the latest code from Github (which, I assume, doesn't contain the changes you're talking about). I used the same solution script as in the OP (just replaced "vertices of group 3" with "vertices of group 2"), and picked points randomly from the quadrilateral with corners at (0, 0), (0, 1), (-0.1, 4), (-1, 4). Using PointsProbe with close_limit=0.0, out of the 92475 random points as many as 55155 returned NaN.

I hope this helps. =)

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.