ufl team mailing list archive

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

Hi,

I have some problems with the built-in function Constant() when generating code
using the quadrature option. The problem is that it is restricted by default to
('+/-'). Could someone remind me why this is the case?

The reason this causes problems is that the map for restriction.CONSTANT
(ufcformat.py) is "0". So if one does:

el = FiniteElement("Lagrange", "triangle", 1)
DG = FiniteElement("DG", "triangle", 0)

v = TestFunction(el)
u = TrialFunction(el)
C = Constant("triangle")
Cu = C*u

# 1
a = dot(grad(v), grad(Cu))*dx

Transforms of the kind Jinv0_00 will appear in the header file. These transforms
are of course not defined for a cell integral which will result in a compile
error in DOLFIN. Tensor representation works because terms involving Jinv0_00
fall away.

A second form that will fail (at FFC compile time) is:

# 2
a = dot(jump(v), jump(Cu))*dS

Because 'C' is already restricted.

However
# 3
a = C*dot(jump(v), jump(u))*dS

will work.

If 'C' is defined as
C = Function(DG)
things work because it is no longer restricted and the map for None is "", so we
only have transforms of the kind Jinv_00.

Now, why would someone take derivatives of a constant (like in #1) since the
following form is equivalent (and works)?

# 4
a = C*dot(grad(v), grad(u))*dx

Well, for very complex forms involving a lot of functions and constants it is
just easier to let FFC worry about which terms vanishes and which don't. An
optimisation of quadrature would naturally be to remove terms where the
basisfunctions are zero, this will take care of the problem as terms involving
Jinv0_00 disappears.

But since form #2 also fails I think that it would make more sense not to
restrict Constant() by default and let the user handle restrictions.

The way things are the user has to define forms like #3 and #4 or define
constants like:

C = Function(DG)

but then, why do we support Constant()?

What do You have in mind for UFL?

Kristian