Plot of projection onto P1+ cubic bubble looks wrong

[Garth Wells]

Hello,

I am facing a strange issue. I am working on the linear space enriched with the bubble space, V=P1+B.
In this short example, here is the problem:
-First, when I project the function, I get a strange result (the height of my square goes from 1 to 1.6 ...)
-Then, given that bubble function, the integration seems to be wrong (it is like I have a projected function in the integral...)

As I'm working with polynomials, of order 3, I expected the assemble function to be exact. Then, another critic question appears: if this integral is false, it will be the same in the solver, when assembling the form?...

Here is my piece of code:

from dolfin import *

#Spaces
mesh=Rectangle(-2,-2,2,2,20,20)
P1=FunctionSpace(mesh,'CG',1)
B=FunctionSpace(mesh,'Bubble',3)
V=P1+B

#Visualisation Space
mesh2=Rectangle(-2,-2,2,2,100,100)
VL=FunctionSpace(mesh2,'CG',1)

class Image1(Expression):
    def eval(self,value,x):
        if ((x[0])<=1) and (abs(x[1])<=1) and (x[0]>=-1):
            value[0]=1
        else:
            value[0]=0

I_1=Image1()

I=interpolate(I_1,V)

H=project(I,V)

print(assemble(H*dx))

H_visualisation=project(H,VL)

print(assemble(H_visualisation*dx))

plot(H)
plot(H_visualisation)
interactive()

Thanks a lot for any explanation.

Emmanuel

[Andrew McRae (andymc)]

I was about to file a bug about the same thing, then I realised what is happening:

Recall that we can express a FE function u = sum(u_i * phi_i), where phi_i are the basis functions. In a normal function space, there is a set {x_i: i = 1, ..., #dof} such that phi_j = 1 at x_j, and 0 at all the other x_i. In FEniCS, interpolating f onto a function space sets the coefficient u_i = f(x_i). Then u(x_i) = f(x_i) for each x_i (though it may differ between the x_i)

For an enriched function space, this "phi_j = 1 at x_j, and 0 at all other x_i" no longer holds, since the basis functions contain basis functions for (in your case) P1, as well as basis functions for the cubic bubble. So interpolating onto V is equivalent to "interpolate onto P1", "interpolate onto B", and adding the two, which corresponds to what you see.

I'll let someone wiser decide whether this is a bug or intended behaviour...