ffc team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting Anders Logg <logg@xxxxxxxxx>:

> 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;

This code is a copy/paste of what is inside FIAT, so if you can come up with a
simple fix to make evaluate_basis do what you want you should also add it to
FIAT. Personally, I'm not sure it makes much sense to evaluate the basis
outside the element, then you could in principle evaluate it outside the
computational domain (where function spaces and their basis cease to exist in
my world) and still you would expect the result to be *valid*?

Kristian

> 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?
>
> --
> Anders
>