Subject: [FFC-dev] "Bug" in evaluate_basis

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 15:59:49 +0200

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

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 16:03:39 +0200

> functions outside the element?

How is this possible? The functions are not defined outside the element.

Garth

> --
> Anders
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev

---

From: Matthew Knepley <[email-masked]>
Date: Wed, 16 Sep 2009 09:05:20 -0500

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 <[email-masked]> wrote:

> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkqw71UACgkQTuwUCDsYZdEIagCeK7i8b18hMN/aD1tQioosHu6L
> GVsAn28SJ86qU+zUm2AUuYLpvTgTD7ds
> =qOFA
> -----END PGP SIGNATURE-----
>

-- 
What most experimenters take for granted before they begin their experiments
is infinitely more interesting than any results to which their experiments
lead.
-- Norbert Wiener

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 16:24:58 +0200

@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

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 16:34:29 +0200

> they are also a basis for P1 on R^2. You can evaluate 1 - x - y for

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.

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 16:49:56 +0200

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

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 17:37:13 +0200

> evaluate_basis. If you do (1, 1 + eps), then you get -1 - eps.

The sparsity is a result of the basis being non-zero only locally, so I 
would call getting zero outside of the cell from ufc::evaluate_basis(..) 
a feature.

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 17:44:00 +0200

> would call getting zero outside of the cell from ufc::evaluate_basis(..)

That's true only for DG elements. All other bases we know of have
larger support (on a patch of elements).

Furthermore, evaluate_basis does currently *not* return zero outside
the cell, only for *certain* points outside the cell.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 17:50:05 +0200

> larger support (on a patch of elements).

We can define the basis on a cell to be zero outside of the cell. The 
dof map takes care of patching together the functions on neighbouring cell.

> the cell, only for *certain* points outside the cell.

It would be nice if it returned zero outside the cell, but it is a 
relatively low-level function so there may be efficiency reasons why 
this isn't desirable, in which case it wouldn't both me if the returned 
values for points outside the cell are meaningless.

Garth

---

From: Matthew Knepley <[email-masked]>
Date: Wed, 16 Sep 2009 10:53:41 -0500

I do agree with Anders, that if we associate a single basis function to each
vertex (P1),
then that function must have support on the "star" of the vertex, or every
cell which
includes that vertex on its boundary.

  Matt

On Wed, Sep 16, 2009 at 10:50 AM, Garth N. Wells <[email-masked]> wrote:

> We can define the basis on a cell to be zero outside of the cell. The dof
> map takes care of patching together the functions on neighbouring cell.
>
> relatively low-level function so there may be efficiency reasons why this
> isn't desirable, in which case it wouldn't both me if the returned values
> for points outside the cell are meaningless.
>
> Garth
>
>  --
>> Anders
>>
>
>
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 17:56:01 +0200

> values for points outside the cell are meaningless.

Yes, that could be an option, or adding to the manual that the values
are undefined outside the cell.

It's just that it would be very useful for us (Marie and I) to reuse
the cell basis functions on cell patches. We need it for the automatic
error estimation (which is close to working btw).

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 17:58:11 +0200

You mean the transitive closure of the support? ;-)

--
Anders

On Wed, Sep 16, 2009 at 10:53:41AM -0500, Matthew Knepley wrote:
> then that function must have support on the "star" of the vertex, or every cell
> which

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 18:02:51 +0200

Anders Logg wrote:
> the cell basis functions on cell patches. 

Are you wanting to evaluate functions like v2 = x outside of [0, 1]? 
What do you do with it? I added a Blueprint to UFC/FFC to evaluate basis 
functions at reference coordinates. Would that help with reuse?

Garth

We need it for the automatic

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 18:07:49 +0200

On Wed, Sep 16, 2009 at 06:02:51PM +0200, Garth N. Wells wrote:
> What do you do with it?

I want to have a basis for say P2 on a patch of elements surrounding an
element K and one option would be to use the extension of the P2
basis on K.

> I added a Blueprint to UFC/FFC to evaluate basis

You mean evaluating at a point not on the physical element but on the
reference element? Like evaluate_basis but without the cell argument?
No, that would not help as those points would also need to avoid the
special line where the values are undefined.

--
Anders

---

From: Matthew Knepley <[email-masked]>
Date: Wed, 16 Sep 2009 11:11:20 -0500

> You mean the transitive closure of the support? ;-)

If you express it as a sieve, Yes.

   Matt

> >                         @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
> >                 The sparsity is not a result of evaluate_basis returning
> 0
> >                 outside the
> >                 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
> >                 For example, if you take the first basis function (1 - x
> - y)
> >                 on the
> >                 expected although that point is outside the triangle. But
> if
> >                 you try
> >             The sparsity is a result of the basis being non-zero only
> locally,
> >             so I
> dof
> this
> iEYEARECAAYFAkqxCxMACgkQTuwUCDsYZdEm6ACeIwB2B5lEezV67AtPczQtdCnd
> btgAn0ylMIUDs49LTyM1yv+W72+ISeXM
> =QSEr

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 18:12:48 +0200

> element K and one option would be to use the extension of the P2

I don't understand how evaluating the P2 basis for K outside of K would 
work. For example, it could be greater than one.

> reference element? Like evaluate_basis but without the cell argument?

Yes.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 16 Sep 2009 18:26:26 +0200

>     x = 2.0 *x/(1.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

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 18:27:26 +0200

> work. For example, it could be greater than one.

Yes, it can be greater than one but that's not a problem. The point is
that the extension of the basis for P2 on K spans P2 on R^2.

Just think of the nodal basis for P1 on the reference triangle:

  v1 = 1 - x - y
  v2 = x
  v3 = y

These three functions span P1 on R^2 or on any other triangle than K.
But they are only a *nodal* basis on K.

--
Anders

---

From: Shawn Walker <[email-masked]>
Date: Wed, 16 Sep 2009 12:57:52 -0400 (EDT)

> But they are only a *nodal* basis on K.

This discussion doesn't concern me, but here is my opinion.

I think that it should be possible for a user to evaluate the basis 
functions outside the reference element.  This could potentially be useful 
for a generalized FEM.  Or perhaps a user wants to `hack' DOLFIN to do 
some peculiar thing, for whatever reason.  This should be allowed!

I think in general, one should not limit functionality unless there is a 
strong reason otherwise.  You could still output a warning message.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 19:07:14 +0200

On Wed, Sep 16, 2009 at 12:57:52PM -0400, Shawn Walker wrote:
> functions outside the reference element.  This could potentially be
> useful for a generalized FEM.  Or perhaps a user wants to `hack'
> DOLFIN to do some peculiar thing, for whatever reason.  This should
> be allowed!
>
> I think in general, one should not limit functionality unless there
> is a strong reason otherwise.  You could still output a warning
> message.

A share this view.

In this particular case, the strong reason for disallowing it might be
that it's difficult to get around because of the particular way FIAT
has implemented the basis functions (by a singular map from a square).

If that's the case, then we can work around it in our code (by
expanding in the Bernstein basis and then restricting back to the
cell).

Rob?

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 16 Sep 2009 19:58:26 +0200

>> is a strong reason otherwise.  You could still output a warning

What is the isoparametric map outside of the cell?

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 16 Sep 2009 20:02:45 +0200

> What is the isoparametric map outside of the cell?

The map is a polynomial and so are the basis functions so they
naturally extend to any point in the universe.

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 22 Sep 2009 17:21:56 +0200

> naturally extend to any point in the universe.

I've added a note to the UFC manual about the value being undefined at
a point outside the cell. It's the easiest solution...

--
Anders