I have a rather simple script to solve the Poisson equation with a source and a sink term on a certain boundary. After calculating the steady-state solution, I would like to compute some derived quantities from it. I understand that for saving these quantities to a file, I need to convert them from their UFL representation by either using the `project`

or `interpolate`

function.

```
from fenics import *
import numpy as np
# Define tolerance
tol = 1E-14
# Create mesh and define function space
mesh = IntervalMesh(100, -0.0005, 0.001)
V = FunctionSpace(mesh, 'P', 1)
# Define subdomains for regions c and a
subdomain_c = CompiledSubDomain('x[0] <= mesh_boundary + tol', mesh_boundary=0., tol=tol)
subdomain_a = CompiledSubDomain('x[0] >= mesh_boundary - tol', mesh_boundary=0., tol=tol)
# Define regions to set coefficients on
regions = MeshFunction('size_t', mesh, mesh.topology().dim())
subdomain_c.mark(regions, 0)
subdomain_a.mark(regions, 1)
# Define coefficient class
class Variable_Coefficient(UserExpression):
def __init__(self, regions, c_0, c_1, **kwargs):
super().__init__(**kwargs)
self.regions = regions
self.c_0 = c_0
self.c_1 = c_1
def eval_cell(self, values, x, cell):
if self.regions[cell.index] == 0:
values[0] = self.c_0
elif self.regions[cell.index] == 1:
values[0] = self.c_1
# Assign coefficient values
D = Variable_Coefficient(regions, 1E-6, 1E-7, degree=0)
G = Variable_Coefficient(regions, 1, 0., degree=0)
R = Variable_Coefficient(regions, 0., 0.5, degree=0)
u_max = Constant(2.5)
u_sink = Constant(1.)
# Define variational problem
u = Function(V)
v = TestFunction(V)
S = conditional(lt(u, u_sink), u, u_sink)
F = D * dot(grad(u), grad(v)) * dx - G * u_max * v * dx + G * u * v * dx + R * S * v * dx
# Compute solution
solve(F == 0, u)
# Compute decomposition and redeposition rates
decomposition = project(G * u_max - G * u, V)
redeposition = project(R * S, V) #this works
#decomposition = interpolate(G * u_max - G * u, V)
#redeposition = interpolate(R * S, V) #this does not work
# Save solution to file in VTK format
vtkfile = File('/mnt/c/Users/margenfeld/Desktop/fenics/output/solution.pvd')
vtkfile << u
vtkfile << decomposition
vtkfile << redeposition
```

This works with `project`

, however, I get some aphysical oscillations in the saved quantities. Thus, I wanted to check if this is also the case for `interpolate`

. When I run the program, I get the error message

```
Traceback (most recent call last):
File "/home/user/fenics/fin_decomposition.py", line 76, in <module>
dec = interpolate(decomposition, V)
File "/home/user/anaconda3/envs/fenics/lib/python3.8/site-packages/dolfin/fem/interpolation.py", line 73, in interpolate
Pv.interpolate(v)
File "/home/user/anaconda3/envs/fenics/lib/python3.8/site-packages/dolfin/function/function.py", line 365, in interpolate
self._cpp_object.interpolate(u)
AttributeError: 'Sum' object has no attribute '_cpp_object'
```

I was of the impression that `project`

and `interpolate`

can be used interchangeably in terms of syntax. That does not seem to be the case, however. I am looking forward to a solution to this problem!