Hi,

I am trying to compute sensitivities using the adjoint method in DOLFIN (without using dolfin-adjoint), similar to this post.

My adjoint equation is as follows:

\frac{\partial{R}}{\partial{U}}\cdot \lambda^T = \frac{\partial{J}}{\partial{U}}, where \frac{\partial{J}}{\partial{U}} = \frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}}

However I am struggling to assemble the RHS of the adjoint equation (\frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}}) in a consistent manner to solve to adjoint equation, resulting in a `TypeError: unsupported operand type(s) for *: 'Form' and 'Form'`

error.

Please see adapted MWE below:

```
import numpy as np
from dolfin import *
from dolfin_adjoint import *
n = 10
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 1)
S = FunctionSpace(mesh, "DG", 0)
u = Function(V) # test function (displacment field)
v = TestFunction(V) # Test function
v_ = TrialFunction(V) # trial function
# kinematics
lambda_ = 0.5769
mu = Function(S)
mu.vector()[:] = 0.3846
P = 2.0*mu*sym(grad(u)) + lambda_ * \
tr(sym(grad(u)))*Identity(len(u))
dx = Measure("dx", metadata={"quadrature_degree": 4})
R = inner(grad(v), P)*dx
# def left(x): return np.isclose(x[0], 0)
# def right(x): return np.isclose(x[0], 1)
def left(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0, tol)
def right(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 1, tol)
bc_L = DirichletBC(V, Constant((0,0)), left)
bc_R = DirichletBC(V, Constant((1,1)), right)
bcs = [bc_L, bc_R]
J = derivative(R, u, v_)
# solve the problem
problem = NonlinearVariationalProblem(R, u, bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
# set an arbitrary functinal
psi = (mu*dot(u, u))
J = (psi**2)*dx
# Derivative of J w.r.t psi
dJ_dpsi = 2*psi*dx
# partial derivative of J w.r.t. mu
dJdmu = assemble(derivative(J, mu))[:]
# partial derivative of R w.r.t. mu
dRdmu = assemble(adjoint(derivative(R, mu))).array() # partial derivative
lmbda = Function(V)
# solve the adjoint vector
lhs = adjoint(derivative(R, u))
# RHS using chainrule
rhs = -(dJ_dpsi * derivative(psi*dx, u)) # TypeError: unsupported operand type(s) for *: 'Form' and 'Form'
problem = LinearVariationalProblem(lhs, rhs,lmbda, bcs=bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
lmbda_vec = lmbda.vector()[:]
dJdmu += np.dot(dRdmu,lmbda_vec)
```

Is there way to perform the \frac{\partial{J}}{\partial{\psi}}\frac{\partial{\psi}}{\partial{U}} calculation such that the output (`rhs`

) can be used in the solver?

I have tried assembling the forms and then performing the multiplication (`assemble(dJ_dpsi) * assemble(derivative(psi*dx, u)`

), however this is not compatible with the solver inputs since the rhs is now a `vector`

rather than `form`

.

Thanks!