Some problems with quadrature expressions

[Nico Schlömer]

When defining and Expression, by default it gets immediately interpolated to a polynomial function of sorts. Since 58f6362, this needs to be made explicit via

Expression('1/x[0]', degree=5)

(or even specifying the element), so nobody can say that they didn't know. So far, so good.

Once in a while, I encounter rather curious expressions like

Expression('(x[0]**2 + x[0] * sin(x[0]) + cos(x[0]) - 1) / (x[0] * x[0])', degree=5)

This looks like it has a pole at 0, but really it doesn't: it all cancels out nicely. Of course, eval doesn't know about that and just spits out NaNs if you try to evaluate it at x[0]==0. To avoid this, it is necessary (see here) to declare the Expression as a "quadrature element":

Q = FiniteElement('Quadrature', triangle, degree=degree, quad_scheme='default')

Unfortunately, you already have to specify the geometry here (triangle), but I guess that's because already now, a quadrature scheme for the element is set in stone. This leads to further complications: You can't just

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
v = TestFunction(V)

Q = FiniteElement("Quadrature", triangle, degree=5, quad_scheme='default')
f = Expression('1.0 / x[0]', element = Q)

b = assemble(f * v * dx)

but you'll have to use the exact quadrature degree from f:

b = assemble(f * v * dx(degree=5))

Meaning: Wherever f is integrated, you'll have to have the corresponding degree available and explicitly specify it. This messes up parts of my code since I have to pass around another variable all the time. Even worse: If there's a second Expression with a different degree, integration is flat out impossible:

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
v = TestFunction(V)

Q = FiniteElement('Quadrature', triangle, degree=5, quad_scheme='default')
f = Expression('sin(x[0]) / x[0]', element=Q)

P = FiniteElement('Quadrature', triangle, degree=3, quad_scheme='default')
g = Expression('sin(x[0]) / x[0]', element=P)

assemble(f * g * v * dx(degree=3))

One possibility how the above could be improved is the following:

When creating a quadrature expression, the desired degree is specified, but the quadrature rule isn't set. When it comes to integrating, one checks the combined degree of the preceding expression (e.g., degree 9 in the case of f * g * v above), and the rule is chosen according the the geometry.

Opinions?

[Jan Blechta]

This is hardly a DOLFIN issue and very probably won't be fixed until the release. @martinal?

[Martin Sandve Alnæs]

My opinion: No.

In fact I deliberately removed that kind of behaviour because the implementation had lots of buggy side effects.

The f and g objects in your code needs to have a fully specified element associated for dolfin assembly to evaluate them in the right points. Thus the degree cannot just be chosen internally in the form compiler code.

To get the behaviour you want today, use UFL expressions (instead of Expression) which gets inlined in the form and computed exactly at quadrature points.

Situations exist where UFL expressions are not sufficient, but you haven't given any reason why they're not sufficient for you. I believe there are long standing issues in bugtrackers discussing alternative future solutions, but so far nobody has cared enough to do the work.

[Nico Schlömer]

Thanks everyone for the feedback.

To get the behaviour you want today, use UFL expressions (instead of Expression) which gets inlined in the form and computed exactly at quadrature points.

Sounds just like what I need. Are there examples out there? (Never seen those in the wild.)

I believe there are long standing issues in bugtrackers

I'd be great to have a link to one of those so I can read up on the discussion.

[Nico Schlömer]

The f and g objects in your code needs to have a fully specified element associated for dolfin assembly to evaluate them in the right points.

Aha. You mean that the polynomial degree alone is not enough to deduct a quadrature rule? (I'm always thinking in terms of Lagrange elements so I hadn't thought this could be an issue.)

[Martin Sandve Alnæs]

No that's not at all what I mean.

From what I understand, you're suggesting creating a dolfin.Expression object (which is a ufl.Coefficient) with a quadrature element without specifying the quadrature rule. It's not possible to evaluate that object in quadrature points when the rule is not specified. Dolfin simply asks that object to restrict itself to an element, so the object must know which points to evaluate itself in. That information must be inherent in the object.

But sorry, I don't have more time for this discussion.

[Martin Sandve Alnæs]

Just set x = SpatialCoordinate(mesh) and use any UFL operators on that.

[Jan Blechta]

I have to pass around another variable all the time

You can fetch the information you need to build the integration measure from the integrand if it already contains a quadrature element coefficient.

Unfortunately, you already have to specify the geometry here (triangle)

Supporting expressions with no cell/gdim specified is rather a historical heritage. It just brings a lot of pain throughout whole FEniCS codebase.

Q = FiniteElement('Quadrature', triangle, degree=5, quad_scheme='default')
f = Expression('sin(x[0]) / x[0]', element=Q)

P = FiniteElement('Quadrature', triangle, degree=3, quad_scheme='default')
g = Expression('sin(x[0]) / x[0]', element=P)

assemble(f * g * v * dx(degree=3))

The reason why this can't work is not purely technical but it is a matter of a definition of quadrature element in my opinion. Quadrature element means "take quadrature of an integrand using a specified quadrature rule". (Other way to look on this is that a function on a quadrature element is a quadrature rule with a weight.) One cannot require to take quadrature of f*g*v using Q and using P at the same time. Giving it any special meaning (like the one you suggested) does not sound consistent.

What you're asking for, I think, are actually pointwise expression (or however it is called). What you want is integrate f*g*v evaluating f and g at (some) quadrature points (without intermediate restriction to element) taking some quadrature rule of degree deg(f) + deg(g) + deg(v) = 5 + 3 + 1 = 9. That's just different functionality than Quadrature element but related. Note that this has been discussed already somewhere but I'm not sure if there is an open blueprint.

On the other hand the following code might work, but it's probably too technically complicated to make it work in the current FFC.

Q = FiniteElement("Quadrature", triangle, degree=5, quad_scheme='default')
f = Expression('1.0 / x[0]', element = Q)

b = assemble(f * v * dx)  # Infer the quadrature rule from f