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]