Marie and I came across a "bug" when trying to use evaluate_basis in a
way it was probably not intended for.
In particular, we want to evaluate the basis functions for a cell in a
point that lies outside the cell. This fails when we evaluate the
basis functions in a point where the collapsing square map in FIAT is
singular.
Here's a simple script that illustrates the problem:
from dolfin import *
import numpy
mesh = UnitSquare(1, 1, "left")
V = FunctionSpace(mesh, "CG", 1)
element = V.dolfin_element()
cell = Cell(mesh, 0)
value = numpy.zeros(1)
# Value should be 1.0 (constant along x = 1)
element.evaluate_basis(1, value, numpy.array((1.0, 1.0000000000001)), cell)
print value
# Should also be 1.0 but result is 0 (eps)
element.evaluate_basis(1, value, numpy.array((1.0, 1.0000000000000)), cell)
print value
The problem is the following lines in evaluate_basis:
// Map coordinates to the reference square
if (std::abs(y - 1.0) < 1e-14)
x = -1.0;
else
x = 2.0 *x/(1.0 - y) - 1.0;
y = 2.0*y - 1.0;
In the above example, y = 1 so x is mapped to -1.0 on Rob's reference
square.
Is there a simple fix that would allow us to evaluate the basis
functions outside the element?
