On Wed, Sep 16, 2009 at 04:34:29PM +0200, Garth N. Wells wrote:
Anders Logg wrote:
@Garth: Yes, they are defined everywhere. For example, v1 = 1 - x - y,
v2 = x and v3 = y are a nodal basis for the reference triangle, but
they are also a basis for P1 on R^2. You can evaluate 1 - x - y for
any values of x and y.
I don't agree at all. The function 1-x-y can be evaluated anywhere, but
the basis functions are defined only on the cell. If the basis functions
were defined everywhere, the method would loose all sparsity.
The sparsity is not a result of evaluate_basis returning 0 outside the
cells, it's a result of the assembly process and the local-to-global
mapping.
Either way, that's not how evaluate_basis is implemented today. It
works perfectly fine to evaluate the natural extensions of the basis
functions, except along one line where the mapping happens to be
singular. This is not in anyway connected to how the basis functions
should be defined, it's just a consequence of the particular way Rob
has implemented the basis functions in FIAT (via a mapping from Jacobi
polynonials on a square).
For example, if you take the first basis function (1 - x - y) on the
reference triangle and evaluate it at (2, 2), you get -2 as
expected although that point is outside the triangle. But if you try
to evaluate it at (1, 1), you get 0 because of a bug/feature in
evaluate_basis. If you do (1, 1 + eps), then you get -1 - eps.