ffc team mailing list archive

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

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

Attachment: signature.asc
Description: Digital signature