GridFunction: boundary dof evaluation for 2D mesh

Dear Bempp community,

I am observing a strange behaviour when using the GridFunction method to project a function onto a 2D mesh.

In the example below, I project the same function onto a 2D square mesh and onto a 3D cube mesh. In theory, the square mesh should coincide with the bottom face of the cube. However, the projection on the 2D square mesh appears to be shifted or incorrect along the boundary, as shown in the attached image.

Does anyone know why this might be happening, and how to avoid or correct this issue?

Thank you in advance — I look forward to your advice.

EDIT: I think it’s to do with P1 function space as DP1 and DP0 were fine. Can we still use P1 for 2D mesh?

import bempp_cl.api as bempp
import numpy as np

x,y,h=1,1,0.1
freq=50
omega = 2*np.pi*freq
c=340
k = omega/c

d = np.array([1., 1., 1])
d /= np.linalg.norm(d)

@bempp.complex_callable
def u_inc(x, n, domain_index, result):
    result[0] = np.exp(1j * k * np.dot(x, d))

grid_2d=bempp.shapes.screen(np.array([[0,0,0],[0,y,0],[x,y,0],[x,0,0]]),h=h)
grid_3d=bempp.shapes.cube(1,[0,0,0],h=h)

p1_space_2d=bempp.function_space(grid_2d,'P',1)
p1_space_3d=bempp.function_space(grid_3d,'P',1)
pd_2d=bempp.GridFunction(p1_space_2d,fun=u_inc)
pd_3d=bempp.GridFunction(p1_space_3d,fun=u_inc)
pd_2d.plot()
pd_3d.plot()