Miklos is exactly right. What you see in the very corners is the following: if you have a quad with 3 zeroes and a nonzero (let's say a 1), if the quad is split along a 0-0 diagonal the you get a 0-0-0 triangle and a 0-1-0 triangle, but if the quad is split a 0-1 diagonal then you get two 0-0-1 triangles.

It's even worse if the quad is
+1 -1
-1 +1

in which case you get a 'valley' or a 'ridge', even though the actual function is symmetric about 0.

To answer your last question, the triangular mesh isn't symmetric in that way.

On 1 August 2016 at 11:08, George Ovchinnikov <lives9@gmail.com> wrote:

Dear all,

I'm solving Poisson equation and I found some irregularity. Here the code:

from firedrake import *

mesh = UnitSquareMesh(8, 8, quadrilateral=True)

V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 0)
w = Function(W)
w.interpolate(Expression(1)) #Initial material distribution

bc0 = DirichletBC(V, Expression(("0")), 1)
bc1 = DirichletBC(V, Expression(("0")), 2)
bc2 = DirichletBC(V, Expression(("0")), 3)
bc3 = DirichletBC(V, Expression(("0")), 4)
bcs = [bc0, bc1, bc2, bc3]

u = TrialFunction(V)
v = TestFunction(V)

f = Function(V)
f.interpolate(Expression("1"))
L = f * v * dx
rhs = assemble(L)

a = (w * dot(grad(u), grad(v))) * dx
frm = assemble(a, bcs = bcs)
sol = Function(V)
solve(frm, sol, rhs, solver_parameters = {'ksp_type': 'preonly',
'pc_type': 'lu'})
File("u.pvd") << sol

The asymmetry of solution (see the figure attached) is clearly visible:
for example, upper left element is not the same as upper right.

I'm wondering what causes this and can this be remedied?

Also, I'm getting exactly same result (with the same asymmetry) on
triangular, not quadrilateral mesh. Why solutions of triangular and
quadrilateral meshes are the same?

Sincerely,
George