Hi Tuomas,
It's possible that you're hitting a non-affine issue here. At the moment, the Jacobian, dx_i/dX_j, of the coordinate transformation from the reference cell to the physical cell, x(X), is treated as being constant over the cell. For non-simplex cells (i.e. quadrilaterals, triangular prisms, tetrahedra), this is only an approximation. The Jacobian matrix is used when replacing the measure (dx == |det J|*dX) and in derivatives (d/dx_i == sum_j dX_j/dx_i d/dX_j).
I have some work in progress, hopefully landing in the next few weeks, after which the Jacobian will be calculated separately at each quadrature point. This is currently spread across branches in different components of Firedrake. I'll try your code tomorrow to see if the results come out the same (or at least much closer!).
Just checking, I assume that doing the problem on tetrahedra gives matching results, even with a deformed mesh? (If so, this suggests further that it's a non-affine approximation issue). I.e., replace the first four lines of code with