 Interpolate() function not working properly with BDM elements

Hello, I am using FEniCS 2019.01 and am experiencing and issue when trying to interpolate a specific Expression to the BDM1 space. The problem is set on a isosceles triangle element with vertices (1,0), (0,1), (-1,0). The considered function is v=(0,x^2). As you can see from the output of the code below, the interpolated function seems to be Iv = (\frac{2x}{3}, -\frac{y}{3}+\frac{1}{9}).

Yet when calculating the interpolated function by hand and in NGSolve with the method provided in https://ngsolve.org/docu/latest/i-tutorials/unit-2.10-dualbasis/dualbasis.html, I get to the interpolant Iv = (\frac{x}{2},-\frac{y}{2}+\frac{1}{3}).

So is there an issue here in FEniCS, or am I making a mistake?

from dolfin import *

editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, 'triangle', 2, 2)
editor.init_vertices(3)
editor.init_cells(1)
editor.close()
mesh.init()
x = SpatialCoordinate(mesh)

W_BDM = FunctionSpace(mesh, "BDM", 1)

v = Expression(('0','x*x'), degree = 2)
Iv_guess = as_vector([2.0*x/3.0,-x/3.0+1.0/9.0])

Iv = interpolate(v, W_BDM)

print('L2-norm of Iv-Iv_guess:', str(sqrt(assemble((Iv-Iv_guess)**2 * dx(mesh)))))

output:

L2-norm of Iv-Iv_guess: 2.127340761518145e-16

After more reflection, I really think there is something wrong with the BDM interpolation using the interpolate() function:

The BDM interpolation operator is supposed to map divergence-free functions to divergence-free BDM functions. As you can see in the above example, the function which FEniCS calculates is clearly not divergence-free even though the input function is.

Is my observation correct and is there a way to fix this?