Hi,

I have reworked the Poisson solver from the FEniCS tutorial, to solve an optimal control problem constrained by a Poisson equation. I am having trouble creating a .vtk file for my function u, which describes the control of the problem. I’ve included the entirety of the code below, since the function u is involved in several steps. I get the following error:

```
RuntimeError Traceback (most recent call last)
/tmp/ipykernel_530/1586109424.py in <module>
77 vtkfile = File('poisson/optimalcontrol.pvd')
---> 78 vtkfile << u
79 vtkfile = File('CondGDplots/optimalstate.pvd')
80 vtkfile << y
~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/io/__init__.py in __lshift__(self, u)
19 self.write(u[0], u[1])
20 else:
---> 21 self.write(u)
22
23
RuntimeError: Unable to cast Python instance to C++ type (compile in debug mode for details)
```

```
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(16, 16)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Constant(0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Init problem data
x0, x1 = MeshCoordinates(mesh)
yd = conditional(x0>0.5,1,-1)
beta = Constant(10)
lmbd = 0.01
# Initialize u
u = Constant(0)
for i in range(10):
#Solve state y_n
y = TrialFunction(V)
v = TestFunction(V)
f = beta*u
a = dot(grad(y), grad(v))*dx
L = f*v*dx
# Compute solution
y = Function(V)
solve(a == L, y, bc)
#vtkfile = File('CondGDplots/state.pvd')
#vtkfile << y
#Solve adjoint state p_n
p = TrialFunction(V)
v = TestFunction(V)
f = y-yd
a = dot(grad(p), grad(v))*dx
F = f*v*dx
# Compute solution
p = Function(V)
solve(a == F, p, bc)
#vtkfile = File('CondGDplots/adjointstate.pvd')
#vtkfile << p
#break
#Choose step direction
df = lmbd*u+beta*p
w = conditional(df >= 0 ,10,-10)
#Calculate step state y(v)
yw = TrialFunction(V)
v = TestFunction(V)
f = beta*w
a = dot(grad(yw), grad(v))*dx
L = f*v*dx
# Compute solution
yw = Function(V)
solve(a == L, yw, bc)
#Choose step length
dyd = y-yd
dys = yw-y
dus = w-u
g1 = assemble((dyd*dys+lmbd*u*dus)*dx)
g2 = assemble(0.5*(dys*dys+lmbd*dus*dus)*dx)
gminarg = -0.5*g1/g2
alpha = min(max(0,gminarg),1)
print(alpha)
if alpha < 1e-8:
break
u = u+alpha*dus
vtkfile = File('poisson/optimalcontrol.pvd')
vtkfile << u
vtkfile = File('CondGDplots/optimalstate.pvd')
vtkfile << y
```

Thank you in advance!