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()