 # Project(x) works, interpolate(x) does not

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 <= mesh_boundary + tol', mesh_boundary=0., tol=tol)
subdomain_a = CompiledSubDomain('x >= 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 = self.c_0
elif self.regions[cell.index] == 1:
values = 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!

Only selected objects like `Function`s, `Expression`s, and `UserExpression`s can be interpolated using `interpolate`. The `project` function uses L^2 projection, which is oscillatory when projecting a discontinuous function onto a continuous space. Some options for monotone projection would be to project to a discontinuous Galerkin space (e.g., `FunctionSpace(mesh,"DG",0)` for element-wise constants), or to use a lumped-mass projection, which can be implemented a few ways, as demonstrated in this post or this tutorial.

1 Like

Great, thank you a lot! That’s important info that I couldn’t find in the tutorial.