I think you’re missing a fundamental point that @dokken is making. The broken Laplace operator (\nabla^2 defined on each element)
\nabla^2_h : {V}^{h,2} \rightarrow {V}_\text{DG}^{h,0}.
Your MWE is rather confusing so I’ll just write out something simpler
from dolfin import *
mesh = UnitSquareMesh(32, 32)
x = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "CG", 2)
VDG0 = FunctionSpace(mesh, "DG", 0)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = sin(pi*x[0])*sin(pi*x[1])*v*dx
bcs = DirichletBC(V, Constant(0.0), "on_boundary")
u = Function(V)
solve(a == L, u, bcs)
laplace_u = project(div(grad(u)), VDG0)
import matplotlib.pyplot as plt
plt.gcf().add_subplot(1, 2, 1)
plot(u)
plt.title("$u_h(x,y)$")
plt.gcf().add_subplot(1, 2, 2)
plot(laplace_u)
plt.title(r"$\Pi_{V_{DG}^{0,h}} (\nabla^2 u_h(x,y))$")
plt.show()
