Hello everyone,

I am solving a Poisson equation, and that works well, but I also want to take some partial derivatives from the variational form. I am encountering multiple errors and would really appreciate some help. I am taking two derivatives one with respect to K and one with respect to G.

- When I use dK = TrialFunction I get

`NotImplementedError: Cannot take length of non-vector expression.`

It seems to work if I use dK = Function(V), I do not understand why? Similarly for G.

- I can assemble the derivative with respect to K but not with respect to G. When I try to assemble with respect to G I get the following error

`ufl.algorithms.check_arities.ArityMismatch: Integrand arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),) differ from form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None)).`

Is that because G is on the boundary (ds)?

Thank you very much for the help in advance!

```
def solve_steady_state_heat(V,mesh,boundary_markers, K, G):
bmesh = dl.BoundaryMesh(mesh, "exterior", True)
bdim = bmesh.topology().dim()
bot_top_mesh = bot_mesh(mesh,boundary_markers,bmesh,bdim)
boundary_mesh = bmesh #for ease going from 2D to actual mesh
submesh_bottom = bot_top_mesh[0]
submesh_top = bot_top_mesh[1]
# boundary measure and normal vector
# ds(1) --> bottom, ds(2) --> top
ds = dl.Measure("ds", domain=mesh, subdomain_data=boundary_markers)
# define Trial and Test Functions
u_trial = dl.TrialFunction(V)
u_test = dl.TestFunction(V)
# initialize input functions
f = dl.Constant(0.0)
u0 = dl.Constant(253.0)
bmesh_markers = dl.MeshFunction('size_t', bmesh, bdim)
bmesh_markers.set_all(0)
# 1 -- bottom
# 2 -- top
for i, facet in enumerate(dl.entities(bmesh, bdim)):
p_meshentity = bmesh.entity_map(bdim)[i]
p_boundarynumber = boundary_markers.array()[p_meshentity]
bmesh_markers.array()[i] = p_boundarynumber
bottom_submesh = dl.SubMesh(bmesh, bmesh_markers, 1)
top_submesh = dl.SubMesh(bmesh, bmesh_markers, 1)
bc = dl.DirichletBC(V, u0, boundary_markers, 2)
def a(K,G,u_trial,u_test,f):
return dl.inner(ufl.exp(K) * dl.grad(u_trial), dl.grad(u_test)) * dl.dx - f * u_test * dl.dx - G*u_test*ds(1)
res_form = steady_state_heat_varf(Kappa,G,u_trial,u_test,f)
# dK = dl.TestFunction(V) <-- when. I use I get the error
# dG = dl.TestFunction(V)
dK = dl.Function(V) #<--- Question 1) Why does this work? I thought we should use TestFunction
dG = dl.Function(V)
theta_kappa_form = dl.derivative( res_form, K,dK)
theta_kappa = dl.assemble(theta_kappa_form)
theta_G_form = dl.derivative( res_form, G,dG)
theta_G = dl.assemble(theta_G_form)
A_form = ufl.lhs(res_form)
b_form = ufl.rhs(res_form)
A, b = dl.assemble_system(A_form, b_form, bcs=bc)
theta = dl.Function(V)
dl.solve(A, theta.vector(), b)
return theta
```