The following code produces a satisfactory Taylor test:

```
from dolfin import *
from dolfin_adjoint import *
import numpy as np
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
VD = VectorFunctionSpace(mesh, "CG", 1)
W = Function(VD)
W.vector()[:] = .005
ALE.move(mesh, W)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx(mesh)
l = v * dx(mesh)
u = Function(V)
dbc = DirichletBC(V, Constant(1.0), "on_boundary")
# works
# solve(a == l, u, dbc)
# works
A = assemble(a)
b = assemble(l)
dbc.apply(A, b)
solve(A, u.vector(), b)
J = assemble(u ** 2 * dx(mesh))
j = ReducedFunctional(J, [Control(W)])
H = Function(VD)
H.vector()[:] = (np.random.rand(VD.dim())*.5 - 1) * 1e-3
taylor_test(j, W, H)
```

Changing the Dirichlet boundary like below, yields a bad Taylor test:

```
dbc = DirichletBC(V, Expression("x[0]", degree=3), "on_boundary")
```

But with these new boundary conditions and imposing the displacements W, H to be zero on the boundary, we recover the convergence in the residuals:

```
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
VD = VectorFunctionSpace(mesh, "CG", 1)
W = Function(VD)
W.vector()[:] = .005
dbcv = DirichletBC(VD, Constant([0,0]), "on_boundary")
dbcv.apply(W.vector())
ALE.move(mesh, W)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx(mesh)
l = v * dx(mesh)
u = Function(V)
dbc = DirichletBC(V, Expression("x[0]", degree=3), "on_boundary")
# works
# solve(a == l, u, dbc)
# works
A = assemble(a)
b = assemble(l)
dbc.apply(A, b)
solve(A, u.vector(), b)
J = assemble(u ** 2 * dx(mesh))
j = ReducedFunctional(J, [Control(W)])
H = Function(VD)
H.vector()[:] = (np.random.rand(VD.dim())*.5 - 1) * 1e-3
dbcv.apply(H.vector())
taylor_test(j, W, H)
```

From here I am tempted to conlcude that the code:

```
A = assemble(a)
b = assemble(l)
dbc.apply(A, b)
solve(A, u.vector(), b)
```

works for shape optimization, as long as:

- the dirichlet conditions don’t depend on the mesh
- or the boundary of the mesh isn’t moved

Could anybody perhaps further explain what’s going on and in case confirm/deny that my observation is correct? Thanks in advance!

**Addendum**

… and now imposing Neumann conditions on part of the boundary, in the last code snippet, yields the residuals: [1.9974521264767453, 1.9355106278012318, 1.895501764619647], which are close to 2 but decreasing, is this just numerical noise?