Hi, in verifying a very simple “shape optimization” code I have written with dolfin-adjoint, I was performing a test with Hessians and I noted I am not obtaining order 3 residuals. I am posting an example. The overloaded block is just a matrix-vector product (so, 0 Hessian), that enables a transfer from functions on a unit interval to a square.

```
from dolfin import *
from dolfin_adjoint import *
import numpy as np
from pyadjoint import Block
from pyadjoint.overloaded_function import overload_function
# Overloading
def backend_displacement(g, M2, VD):
gv = g.vector()[:]
W = Function(VD)
W.vector()[:] = M2 @ gv
return W
class DisplacementBlock(Block):
def __init__(self, g, M2, VD, **kwargs):
super(DisplacementBlock, self).__init__()
self.kwargs = kwargs
self.add_dependency(g)
self.M2 = M2 # matrix bringing us from g.function_space() to VD
self.VD = VD # the volumetric function space, vector fields version
self.VD_dim = VD.dim() # the volumetric function space, vector fields version
self.V_sph_dim = g.function_space().dim()
def __str__(self):
return 'DisplacementBlock'
def recompute_component(self, inputs, block_variable, idx, prepared):
return backend_displacement(inputs[0], self.M2, self.VD)
def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
adj_input = adj_inputs[0] # this is a vector, with which we test the adjoint jacobian
adj_action = self.M2.T @ adj_input
output = Vector(MPI.comm_world, self.V_sph_dim)
output[:] = adj_action
return output
def evaluate_tlm_component(self, inputs, tlm_inputs, block_variable, idx, prepared=None):
# see http://www.dolfin-adjoint.org/en/latest/documentation/pyadjoint_docs.html
return self.recompute_component(tlm_inputs, block_variable, idx, prepared)
def evaluate_hessian_component(self, inputs, hessian_inputs, adj_inputs, block_variable, idx,
relevant_dependencies, prepared=None):
output = Vector(MPI.comm_world, self.V_sph_dim)
output[:] = 0
return output
radial_displacement = overload_function(backend_displacement, DisplacementBlock)
domain = UnitSquareMesh(5,5)
interval = UnitIntervalMesh(10)
V_int = FunctionSpace(interval, "CG", 1)
V_def = VectorFunctionSpace(domain, "CG", 1)
M2 = np.random.rand(V_def.dim(), V_int.dim())
q = Function(V_int) # null displacement (but doesn't need to be so)
W = radial_displacement(q, M2, V_def)
ALE.move(domain, W)
# PDE
V = FunctionSpace(domain, "CG", 1)
u = Function(V)
w = TrialFunction(V)
v = TestFunction(V)
bcs = DirichletBC(V, 0, "on_boundary")
dx = Measure("dx", domain=domain)
solve(inner(grad(w), grad(v)) * dx == v * dx, u, bcs)
# Cost functional
J = assemble(sin(u) * dx)
Jhat = ReducedFunctional(J, [Control(q)])
d = Function(V_int)
d.vector()[:] = (np.random.rand(len(d.vector())) - .5) * 2 * 1e-2
td = taylor_to_dict(Jhat, q, d)
print(td['R2']['Rate'])
print(td['R2'])
```

The result for the Hessian residual:

```
{'Residual': [6.073355259626869e-09, 1.5186571271107686e-09, 3.7970434087315546e-10, 9.493119657780799e-11], 'Rate': [1.9996975778376476, 1.9998477865513609, 1.9999223191972153]}
```

which is not of third order as expected. The deformation that is being performed only minimally perturbs the mesh, so, no weird compenetrations are happening. What could be wrong then? Thanks a lot in advance