The documentation says the = operator does interpolation. But not clear on what happens when interpolating discontinuous fields. For example:
mesh Th = square(1,1); fespace Vh0(Th,P0); fespace Vh1(Th,P1); Vh0 F0 = x; Vh1 F1 = F0; Vh1 F1p, q; plot(F0,fill=1,value=1,dim=3,wait=1); plot(F1,fill=1,value=1,dim=3,wait=1); solve proj(F1p,q) = int2d(Th)(F1p*q)-int2d(Th)(F0*q); plot(F1p,fill=1,value=1,dim=3);
One way to go from P0 to P1 field is using a projection which gives F1p.
But for doing interpolation (F1) the value should be computed on points. I’m wondering how the “=” operator computes the values at diagonals? Here we have two different values and it doesn’t average them. Is it picking randomly one of the values of F0 to assign F1?