← Back to team overview

ffc team mailing list archive

Re: "Bug" in evaluate_basis

 

@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.

@Matt: We use the local basis on a triangle K as a basis for P2 (and
other spaces) on a patch of elements surrounding K. We make an Ansatz
in terms of the local basis functions. So it's not for extrapolation
of a computed solution.

We could also use the Bernstein polynomials on K or any other
reasonable well-conditioned basis, but it would be convenient to use
the basis that's already in there.

--
Anders


On Wed, Sep 16, 2009 at 09:05:20AM -0500, Matthew Knepley wrote:
> Isn't this strange since the error bounds no longer apply outside the element
> (it is no longer interpolation, but extrapolation)?
>
>   Matt
>
> On Wed, Sep 16, 2009 at 8:59 AM, Anders Logg <logg@xxxxxxxxx> wrote:
>
>     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?
>
>
>     -----BEGIN PGP SIGNATURE-----
>     Version: GnuPG v1.4.9 (GNU/Linux)
>
>     iEYEARECAAYFAkqw71UACgkQTuwUCDsYZdEIagCeK7i8b18hMN/aD1tQioosHu6L
>     GVsAn28SJ86qU+zUm2AUuYLpvTgTD7ds
>     =qOFA
>     -----END PGP SIGNATURE-----
>
>     _______________________________________________
>     FFC-dev mailing list
>     FFC-dev@xxxxxxxxxx
>     http://www.fenics.org/mailman/listinfo/ffc-dev
>
>
>
>
>

> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/ffc-dev

Attachment: signature.asc
Description: Digital signature


Follow ups

References