[firedrake] weirdness on bendy branch

[Stephan Kramer]

Hey guys

More weirdness related to facet integrals from IP viscosity that only 
appears with COFFEE optimisations. This time both on the bendy branches 
and on master (of last week). In the code below, the matrices M, Ms and 
M1+M2 should all be the same.

For DG1 they are not equal on a recent (last week, before any sprint 
merges) install of firedrake/ffc/ufl/coffee on master. It is also 
incorrect on an older install I had, so it's not a recent thing. DG1 
does however produce the right answer on the bendy branches: firedrake 
(bendy_changes) + ffc/ufl (fd_bendy) + coffee (master).

For RT1 on the other hand it's the opposite: the matrices are correct 
and equal on master, but with bendy some of them are different from the 
correct answer.

All of this only if use COFFEE optimisation, i.e. with 
parameters['coffee']={} the matrices are correct and all equal.

I realize at this point the best thing to do is probably to wait if this 
all gets magically fixed when the sprint is over - will report back if 
not - but just thought to let you guys know already

Cheers
Stephan

from firedrake import *
mesh2d = UnitSquareMesh(1,1)
U = VectorFunctionSpace(mesh2d, "DG", 1)
v = TestFunction(U)
u = Function(U)

#parameters['coffee']={}
n = FacetNormal(mesh2d)

def outer_jump(v, n):
     return outer(v('+'), n('+'))+outer(v('-'), n('-'))

F1 =  inner(outer_jump(v, n), outer_jump(u,n))*dS
F2 =  inner(outer_jump(v, n), outer_jump(n,u))*dS
F = inner(outer_jump(v, n), outer_jump(u,n)+outer_jump(n,u))*dS

M1 = assemble(derivative(F1, u)).M.values
M2 = assemble(derivative(F2, u)).M.values
Ms = assemble(derivative(F1+F2, u)).M.values
M = assemble(derivative(F, u)).M.values
print abs(M-(M1+M2)).max()
print abs(M-Ms).max()