Hi,

I’ve got a quick question regarding dolfin-adjoint, in particular the computation of shape derivatives. Assuming that I try to solve the optimization problem of the form

\min J(\Omega) = \int_\Omega (y - y_d)^2 dx

where y_d is some desired state and y(\Omega) is the solution of a PDE on \Omega.

The shape derivative for this is then given by

dJ(\Omega)[V] = \int_\Omega (y - y_d)^2 \nabla \cdot (V) - 2*(y - y_d)* \nabla y_d [V] dx

As far as I understand, if y_d is an ufl expression only depending on a SpatialCoordinate, then the shape derivative will be correct. However, if y_d is a Function or Expression, then the above code will not include the second term.

Is there any way to treat such a problem currectly when using a Function / Expression?

In the following I created a MWE highlighting the problem

```
from fenics import *
import numpy as np
mesh = UnitSquareMesh(25, 25)
dx = Measure('dx', mesh)
x = SpatialCoordinate(mesh)
space = FunctionSpace(mesh, 'CG', 1)
expr = Expression('sin(pi*x[0])*sin(pi*x[1])', degree=2)
y = interpolate(expr, space)
def_space = VectorFunctionSpace(mesh, 'CG', 1)
V = TestFunction(def_space)
y_d_coor = pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)
y_d_expr = Expression('pow(x[0], 2) + pow(x[1], 2) - pow(0.5, 2)', degree=2, domain=mesh)
y_d_func = interpolate(y_d_expr, space)
y_d_list = [y_d_coor, y_d_expr, y_d_func]
J_list = [pow(y - y_d, 2)*dx for y_d in y_d_list]
dJ_ex_list = [pow(y - y_d, 2)*div(V)*dx - 2*(y - y_d)*inner(grad(y_d), V)*dx for y_d in y_d_list]
dJ_mod_list = [pow(y - y_d, 2)*div(V)*dx for y_d in y_d_list]
dJ_ad_list = []
for idx in range(len(J_list)):
vec_ex = assemble(dJ_ex_list[idx])[:]
vec_mod = assemble(dJ_mod_list[idx])[:]
vec_ad = assemble(derivative(J_list[idx], x, V))[:]
error_ex = np.max(np.abs(vec_ex - vec_ad)) / np.max(np.abs(vec_ex)) * 100
error_mod = np.max(np.abs(vec_mod - vec_ad)) / np.max(np.abs(vec_mod)) * 100
print('Relative error: ' + str(error_ex) + ' % (exact)')
print('Relative error: ' + str(error_mod) + ' % (modified)')
print('')
```