[firedrake] Possible bug with perturbed extruded mesh?

[Tuomas Karna]

Hi all,

Encountered this issue with extruded mesh where the coordinates have 
been deformed. The first form a==L just evaluates div(u), the second 
a==L2 is the same but div(u) integrated in parts. Thus the results 
should be equivalent, but that's not what I see:

a == L: 0.999999886488 1.0000001092
a == L2: -0.235329928101 3.47882106918

Without the mesh deformation the two forms give the correct solution 
(1.0). Am I missing something?

Cheers,

Tuomas

---

from firedrake import *
mesh2d = UnitSquareMesh(10,10)
layers = 6
mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers)

P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1)
P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=1)
P1v = VectorFunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1)

# deform mesh
scalar = Function(P1).interpolate(Expression('1.0+x[0]'))
coords = mesh.coordinates
coords.dat.data[:,2] *= scalar.dat.data[:]

u = Function(P1v)
w = Function(P2)
u.interpolate(Expression(('x[0]', '0.0', '0.0')))

tri = TrialFunction(P2)
test = TestFunction(P2)
normal = FacetNormal(mesh)

a = test*tri*dx
L = div(u)*test*dx
L2 = -dot(u, grad(test))*dx + dot(u, normal)*test*(ds_v+ds_t+ds_b)

solve(a == L, w)
print 'a == L:', w.dat.data.min(), w.dat.data.max()

solve(a == L2, w)
print 'a == L2:', w.dat.data.min(), w.dat.data.max()