Hi Jonatan,
it is definitely a bug. As a matter of coincidence, I stumbled upon it two days ago, when creating the new example [1], but did not get to investigating it properly (just worked around the bug poorly by setting close_limit parameter of probes to a higher value, see the example). So thank you for the demonstration scripts! More to the bug that below.
[1] http://sfepy.org/doc-devel/examples/diffusion/time_poisson_interactive.html
On 03/19/2015 06:13 PM, Jonatan Lehtonen wrote:
Hi,
I'm using SfePy to solve a simple Laplace problem in a unit square with Dirichlet boundary conditions (well, what I'm actually doing is a bit more involved, but this is just a minimal example). I'm having some problems with the fact that probing with a PointsProbe occasionally returns nan instead of a proper value. I have attached a mesh file which I crea,ed by using gmsh with the attached .geo file and then manually converting it to the desired format for SfePy (that is, I removed the edges and the z-coordinates of the nodes). The following code, with the attached mesh file, should be able to reproduce this issue:
<snip> > > The code above solves the Laplace problem and creates a PointsProbe with > 100000 random points. It then prints True if any of the probed values equal > nan, and prints the coordinates of the corresponding points. When I run > this code, it usually turns up around 5-10 points which produce a nan > value. This is an example of the output on my system: > <snip>
After some further examination, I've found that the problem occurs in some small triangular regions of the domain. To illustrate this problem, I made a (rather ugly, but possibly helpful) scatter plot of the values in a grid around the first point of the list above, using this (extremely messy) code:
temp = probe.points[np.nonzero(np.isnan(res))[0][0]] X, Y = np.mgrid[temp[0]-0.003:temp[0]+0.003:40j, temp[1]-0.003:temp[1]+0.003:40j] gridprobe = PointsProbe(np.array(np.array([X.flatten(), Y.flatten()]).transpose(), order='C')) plot_data = gridprobe.probe(pb.get_variables()['u'])[1].flatten()
from pylab import * figure() scatter(X.flatten()[~np.isnan(plot_data)], Y.flatten()[~np.isnan(plot_data)]) closest_points = vec.variables['u'].field.get_coor(np.argsort(np.linalg.norm(vec.variables['u'].field.get_coor()
- temp, 2, axis=1))[:4]) scatter(closest_points[:,0], closest_points[:,1], color='r') scatter([temp[0]], [temp[1]], color='g') show()
I have added this:
from sfepy.postprocess.plot_cmesh import plot_wireframe
ax = gca() plot_wireframe(ax, domain.cmesh) ax.axis('equal')
The resulting plot is attached. The triangular area inside the blue grid is where probing results in nans. The green point ('temp' above) is the first nan point which was found earlier (i.e. the point (0.48132644 0.32787828)), and the red points are the four mesh nodes closest to the green point. I included the mesh points because my first assumption was that the nan area would correspond to a mesh triangle, but that does not seem to be the case.
The figure resulting from the above code, zoomed to the region of interest, is attached - it displays also a wire-frame of the mesh.
Now the obvious question at this point is: am I doing something wrong, or is this a bug?
The bug is caused by a wrong assumption I made in past. When evaluating a field in general points, the algorithm, in essence, is as follows: the reference-physical element mapping to the point
- for each point, find the mesh vertex closest to the point's coordinates
- loop over the cells sharing that vertex <--- problem here!!!
- determine the reference element coordinates by applying the inverse of
- if the reference coordinates are within [-close_limit, 1+close_limit], the point is in that cell -> evaluate the FE basis there, and interpolate the field to get the value in the point - otherwise continue 3. if the above fails, the value is undetermined (-> nan in the probe)
As you can see in the attached figure, the points in the white triangle are closest to a mesh vertex, that is shared by four cells (elements) and the cell the point is in is not among them. The points close to the cell edge are blue, because they fall within the close_limit range. Try: gridprobe.set_options(close_limit=0.0).
Regards, Jonatan
P.S. I was originally using evaluate_at() instead of a PointProbe, and that function returns 0.0 instead of nan. I'm not entirely sure if this is intended, but at least I found it extremely confusing.
The probes use evaluate_at() internally, but they set the values in not-found points explicitly to nan, to indicate probing outside of the domain, which can legitimately happen. While evaluate_at() per se was meant to be used inside the domain and should not fail - bad guess :) You are right it is confusing, evaluate_at() should return nan as well.
Until it is fixed properly, which might take some time as the bug is in a low-level cython code, you can hack-around with setting close_limit to a higher value. (Also, it is just the time of grant project proposals and tons of bureaucracy :].)
Could you create a new issue for this, with a link to this discussion?
Thank you for your debugging effort! r.
