 Error creating a vizualizaton .vtk of a UFL function

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, u)
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
L = f*v*dx
# Compute solution
y = Function(V)
solve(a == L, y, bc)
#vtkfile = File('CondGDplots/state.pvd')
#vtkfile << y

p = TrialFunction(V)
v = TestFunction(V)
f = y-yd
F = f*v*dx
# Compute solution
p = Function(V)
solve(a == F, p, bc)
#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
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

Please note that you need to project u and y into a suitable function space for visualization.
Also note that to visualize functions from other spaces than “CG” 1, it is recommended to use XDMFFile.write_checkpoint, as it makes it possible to visualize arbitrary order CG and DG spaces.

1 Like

Thank you for the help! In the initial tests, I did not need to project y before visualization. Which operations did I apply to u that hindered this?

You can only save functions to file.
When you add values to your

You change the type of u.

I appreciate you bear with me. Once again, thanks!