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.