Hi Gianluca,
The solution to the Poisson equation \div \grad u = -1 you show does not seem to obey the boundary condition u = 0.
Compare to:
\[CapitalOmega] =
Region[RegionDifference[Rectangle[{-1, -1}, {1, 1}],
Rectangle[{0, 0}, {1, 1}]]]; Subscript[\[CapitalGamma], D] =
DirichletCondition[u[x, y] == 0, True];uval = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == -1,
Subscript[\[CapitalGamma], D]},
u, {x, y} \[Element] \[CapitalOmega]]; ContourPlot[uval[x, y], {x, y} \[Element] \[CapitalOmega]
