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

I/O
1

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!

2

Only selected objects like Functions, Expressions, and UserExpressions 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.

2 Likes
3

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