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 created 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:

import numpy as np
from sfepy.discrete import (FieldVariable, Integral, Equation, Equations, Problem)
from sfepy.discrete.fem import (Mesh, FEDomain, Field)
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.postprocess.viewer import Viewer


mesh = Mesh.from_file('rectangle2D.mesh')
domain = FEDomain('domain', mesh)
omega = domain.create_region('Omega', 'all')
# Groups 1 and 3 correspond to the bottom and top of the square, respectively
gamma0 = domain.create_region('Gamma 0', 'vertices of group 1', 'facet')
gamma1 = domain.create_region('Gamma 1', 'vertices of group 3', 'facet')
field = Field.from_args('fu', np.float64, 'scalar', omega, approx_order=2)
u = FieldVariable('u', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')
integral = Integral('i', order=2)
term = Term.new('dw_laplace(v, u)', integral, omega, v=v, u=u)
eq = Equation('Laplace equation', term)
eqs = Equations([eq])
ebc0 = EssentialBC('Dirichlet 0', gamma0, {'u.all': 0.0})
ebc1 = EssentialBC('Dirichlet 1', gamma1, {'u.all': 1.0})
pb = Problem('Laplace problem', equations=eqs)
pb.time_update(ebcs=Conditions([ebc0, ebc1]))
vec = pb.solve()
pb.save_state('rectangle2D.vtk', vec)
view = Viewer('rectangle2D.vtk')
view()


from sfepy.discrete.probes import PointsProbe
probe = PointsProbe(np.random.rand(100000,2))
pars, res = probe.probe(pb.get_variables()['u'])
print(np.any(np.isnan(res)))
print(probe.points[np.nonzero(np.isnan(res))[0]])

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:

sfepy: reading mesh [line2, tri3, quad4, tetra4, hexa8] (rectangle2D.mesh)...
sfepy: ...done in 0.03 s
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (13343, 13343)
sfepy: assembling matrix graph...
sfepy: ...done in 0.01 s
sfepy: matrix structural nonzeros: 151301 (8.50e-04% fill)
sfepy: updating materials...
sfepy: ...done in 0.00 s
sfepy: nls: iter: 0, residual: 1.819854e+01 (rel: 1.000000e+00)
sfepy:   rezidual:    0.01 [s]
sfepy:      solve:    0.07 [s]
sfepy:     matrix:    0.01 [s]
sfepy: nls: iter: 1, residual: 4.729175e-14 (rel: 2.598656e-15)
sfepy: point scalars u at [-0.5 -0.5  0. ]
sfepy: range: 0.00e+00 1.00e+00 l2 norm range: 0.00e+00 1.00e+00
sfepy: iconn: 0.016510 s
sfepy: kdtree: 0.000420 s
sfepy: evaluating in 100000 points...
sfepy:   reference field: 3436 vertices
sfepy:   kdtree query: 0.291867 s
sfepy:   ref. coordinates: 0.466726 s
sfepy:   interpolation: 0.135751 s
sfepy: ...done
True
[[ 0.48132644  0.32787828]
 [ 0.1449193   0.69671521]
 [ 0.48135807  0.32673518]
 [ 0.77758596  0.35547628]
 [ 0.64769394  0.79028833]
 [ 0.66355676  0.78836725]
 [ 0.4799345   0.32667195]
 [ 0.48147468  0.3281688 ]]

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()

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.

Now the obvious question at this point is: am I doing something wrong, or is this a bug?

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.