FEM Meshing artifact at nodes with fewer neighbors

I wrote a 2D-FEM solver to solve some diffusion process and wanted to verify my code with a test problem.
The input was $f(x,y) = x^2+y^2$ and I applied the stiffness matrix on it to get $Delta f = 4$.
The lumped mass matrix $M$ and stiffness matrix $D$ are assembled from an analytical solution of triangles.

Unfortunately, I get too high values at nodes with fewer neighbor elements (5 instead of six) than average (see below), whereas elements with more neighbors underestimate the analytical value.
I also tried a finer resolution, but this didn’t help.

<img src="https://i.stack.imgur.com/qsQM4.png" alt="solution of $Delta f = 4$ for $f(x,y) = x^2+y^2$ ” />

Note that I manually added the value 4 at the boundary nodes, since the corresponding entries in the stiffness matrix are 0.

Do you have an idea, what the problem might be?
In my test simulations for the diffusion process, I couldn’t find this meshing artifact. Only in the test problem.