Extract coefficients from forms

Question

Dear FEniCS community,

I would like to extract coefficients (not necessarily constants, possibly expressions) from bilinear/linear forms in order to apply a preprocessing to them before assembling the forms.

For instance, consider the following code

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 2)

expr1 = Expression("x[0]", element=V.ufl_element())
expr2 = Expression("x[1]", element=V.ufl_element())
expr3 = Expression("x[0]", element=V.ufl_element())
expr4 = Expression("x[1]", element=V.ufl_element())

u = TrialFunction(V)
v = TestFunction(V)
a = expr1*inner(grad(u), grad(v))*dx + expr2*u.dx(0)*v*dx + expr3*u*v*dx
f = expr4*v*dx

Given a and f, my final goal would be to extract expr1, expr2, expr3 and expr4, apply some preprocessing on them (obtaining expr1_new, expr2_new, expr3_new and expr4_new) and then assemble the following modified forms:

a_new = expr1_new*inner(grad(u), grad(v))*dx + expr2_new*u.dx(0)*v*dx + expr3_new*u*v*dx
f_new = expr4_new*v*dx
assemble(a_new)
assemble(f_new)

Q1. Is it possible to extract the coefficients and then plug them back in the forms? If so, how? I have had a look at

for integral in a.integrals():
    print integral.integrand()

but I am not sure on how to parse it
Q2. Will this work also for VectorFunctionSpaces, where expressions may be matrices or vectors?
Q3. Do I have a smart way to detect if two or more of the extracted coefficients are the same? For instance, in the example above, I would need to apply my preprocessing only to expr1 and expr2, since the remaining ones are the same.
Q4. Would integration over dx vs ds make any difference in the procedure?

I know that in my example I can just apply my preprocessing to each expression and then assemble directly a_new and f_new. However, the actual need is to apply it e.g. to a form obtained with derivative(), so in that case I would be not writing down the bilinear form explicitly and thus I would really need to get each coefficient from the form.

Thanks,
Francesco

Answer 1

Hello, coefficients method of Form allows you to get the form's coefficients, while replace function of UFL allows you to perform substitution. Consider

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

# Suppose these will be coefficients of your form
coefs = (Expression('1', cell=mesh.ufl_cell(), degree=1),
         Expression(('x[0]', 'x[1]'), degree=1, cell=mesh.ufl_cell()),
         Constant(((1, 2), (3, 4))))

c0, c1, c2 = coefs
form = c0*inner(u, v)*dx+sqrt(c1[0]**2+c1[1]**2)*inner(u, v)*ds + det(c2)*inner(grad(u), grad(v))*dx
# All coefs in the form
print len(form.coefficients())

# Update the coefs: each coef is mapped to `ones` of the corresponding shape
update = dict((c, Constant(np.ones(c.ufl_shape))) for c in form.coefficients())
form_new = replace(form, update)
# All coefs in the form
print len(form_new.coefficients())

# Form with repeated coefficients is identfied correctly
form = c0*inner(u, v)*dx + c0*inner(u, v)*ds
assert len(form.coefficients()) == 1

The example shows that answers to Q2 and Q4 are respectively yes and no. As for Q3, coeffcients can be compared with ==.

Comments

Thanks! I will mark this as correct answer. The last line in your snippet answers the most common Q3 case that I had in mind. Just to clarify for future reference, however, actually expr1 is not == expr3, because they are two different expression (even though they both represent "x[0]" in my example). This limiting case will never happen when using derivative(), so I am comfortable that Q3 has positive answer in all non trivial (non-user provided forms) cases.

Thanks,
Francesco

---

Hi,
I have experimented with your code in some cases, and I would like to ask an additional question.

Suppose that
a = expr1*expr2*inner(grad(u), grad(v))*dx
(or any similar operations between coefficients)

Is there a way to extract the product expr1*expr2 while looping on coefficients, rather than expr1 and expr2 separately?

Thanks

---

I think this is possible with ufl.algorithms in particular iter_expressions and traverse_terminals. Your example expr1*expr2 is a Product instance where none of the operands is an Argument type.

---

Thanks, I managed to put together a working example at
https://gist.github.com/francesco-ballarin/d3198c067fe8193200d85001f20e2ec4

However, this example does not work if I replace

  a = expr3*expr2*(1 + expr1*expr2)*inner(grad(u), grad(v))*dx

e.g. with

 a = expr2*(1 + expr1*expr2)*inner(grad(u), grad(v))*expr3*dx 

(note that in the first added I have moved the multiplication by expr3 to the end). The basic idea of the code is to find a common node of the tree where there are only expressions and no arguments, however it seems that in the latter expr3 is always separated from expr2(1 + expr1expr2)*inner(grad(u), grad(v)).
Are you aware of a better way to trasverse the tree? It is not a big deal in this case, however the replace_in_form function does not seem to work in vector problems, and I suspect that it is related to this, e.g.

V = VectorFunctionSpace(mesh, "Lagrange", 2)

expr1 = Expression("x[0]", cell=mesh.ufl_cell())
expr2 = Expression(("x[0]", "x[1]"), cell=mesh.ufl_cell())
expr3 = Constant(((1, 2), (3, 4)))

u = TrialFunction(V)
v = TestFunction(V)
a = inner(expr3*grad(u), grad(v))*dx + inner(grad(u)*expr2, v)*dx + expr1*inner(u, v)*dx

Thanks!