UFL quadrature degree error

variational formulation
1

Hello,

I’ve got a quick question regarding UFL and the choice of quadrature degrees. Particularly, consider the following problem, which raises an error

from fenics import *

mesh = UnitSquareMesh(10, 10)
dx = Measure('dx', mesh)

assemble(1*dx(degree=3) + 1*dx(degree=4))

due to the different degrees of the integrals.

The background is that I want to do automatic manipulation of form objects, and I could boil it down to the above problem (at least that’s what I think). In particular, my best guess is that the form compiler estimates the degree of the objects automatically which then leads to a situation as above.
Is this intended behavior?

2

You can see a short discussion about this at: https://github.com/FEniCS/ufl/pull/27.

In the case above it does not really make sense to have different integral degree’s. Is the real problem forr assembling a matrix or vector?

3

Thanks for the link, now I understand the problem.

The above code was only some toy example, for the real problems this concerns (mostly) vectors, but also matrices.

As I said I generate the forms automatically (via derivative and replace) and what I’ve found to work is to do the degree estimation before the assembly. So if “F” is the UFL form, I could run

from ufl.algorithms.compute_form_data import  estimate_total_polynomial_degree
assemble(F, form_compiler_parameters= {'quadrature_degree' : estimate_total_polynomial_degree(F)})

which should lead to the “expected” behavior of only having one assemble loop which possibly over integrates some terms, right?

4

That seems right. Does this degree become the largest or smallest degree (in the case of your scalar example)?

5

Well, for the scalar problem it unfortunately does not work, it seems that the degrees given in the measures override the degree specified in assemble

from fenics import *
from ufl.algorithms.compute_form_data import  estimate_total_polynomial_degree

mesh = UnitSquareMesh(10, 10)
dx = Measure('dx', mesh)

form = 1*dx(degree=3) + 1*dx(degree=4)
assemble(form, form_compiler_parameters= {'quadrature_degree' : estimate_total_polynomial_degree(form)})

which still throws the same error.

If I generate the code automatically, I don’t hard code the degree in the measure and in this case I get the maximum as wanted.

6

In general, where will you want to mix quadrature degrees? You can simply split your assemblies into pieces (which in turn can cause a slow down due to multiple loops through cells and kernels).

7

I actually do not want to mix quadratures, but fenics does that somehow, internally.

I just created a MWE for this

from fenics import *
from ufl.algorithms import expand_derivatives
from ufl import replace
from ufl.algorithms.compute_form_data import  estimate_total_polynomial_degree
import numpy as np



mesh = UnitSquareMesh(10, 10)
dx = Measure('dx', mesh)

W = FunctionSpace(mesh, 'CG', 1)
W_vec = VectorFunctionSpace(mesh, 'CG', 1)
V = TestFunction(W_vec)
x = SpatialCoordinate(mesh)

f = Expression('pow(x[0],2) + pow(x[1],2) - 1', degree=2, domain=mesh)
idx = f.id()

u_expr = Expression('sin(2*pi*x[0])*sin(2*pi*x[1])', degree=1)
u = interpolate(u_expr, W)
v_expr = Expression('x[0]*x[0]+x[0]*x[1]+x[1]*x[1]+sin(x[0])+cos(x[0])', degree=1)
v = interpolate(v_expr, W)

F = inner(u,v)*dx - inner(f,v)*dx
dF_ex = inner(u,v)*div(V)*dx - div(f*V)*v*dx
sd_ex = assemble(dF_ex)
dF_ad = derivative(F, x, V)

### First possibility
material_derivative1 = derivative(F, f, dot(grad(f), V))
material_derivative1 = expand_derivatives(material_derivative1)

### Alternative
f_elem = f.ufl_element()
func_space = FunctionSpace(mesh, f.ufl_element())
argument = Function(func_space)
temp = derivative(F, f, argument)
material_derivative2 = replace(temp, {argument : dot(grad(f),V)})

dF_ad1 = dF_ad + material_derivative1
dF_ad2 = dF_ad + material_derivative2

### Fails
# sd1 = assemble(dF_ad1)
# sd2 = assemble(dF_ad2)

### Works
sd1 = assemble(dF_ad1, form_compiler_parameters={'quadrature_degree' : estimate_total_polynomial_degree(dF_ad1)})
sd2 = assemble(dF_ad2, form_compiler_parameters={'quadrature_degree' : estimate_total_polynomial_degree(dF_ad2)})

assert np.allclose(sd_ex[:], sd1[:])
assert np.allclose(sd_ex[:], sd2[:])
assert np.allclose(sd1[:], sd2[:])

The goal is to compute shape derivatives, where u and v are state variables and f should experience a pull-back (similarly to our previous discussion: Shape Derivatives in UFL)

Also, is there any difference between using the replace approach and differentiating directly?

8

Oh, and as a remark: That only happens when I use an expression, if I interpolate the expression to a FE space it works fine. Also using other degrees in the expression totally works fine.

9

I tend to avoid using the Expression class (rather use spatial coordinates, at least for the example you present here), as Expressions are difficult to differentiate correctly.

10

Okay, I see, but still I want to support them.

Thanks a lot for the help.