@TeM I saw the link does not seem to be from the native stack exchange site for image links, if the procedure you take to make the link is tedious note that you can create a link for an image by "pretending" to start to write an answer, copy pasting the picture into your answer and then copy pasting the link that was automatically generated. I found that trick searching on google.

@TeM You should use the same mesh in Plot3D and NDSolve!

@TeM, have a look in the ToElementMesh ref page under the options section.

@TeM Did you read my comment " same mesh in Plot3D and NDSolve". Show your code if you need further assistance!

@UlrichNeumann: Yes and thank you as always. But my ignorance on the subject doesn't allow me to understand how to get out of it. I edited the question by adding more code. Thanks again!

@TeM Your code gives a smooth plot in MMA v12.2

@user21 and Ulrich, this is surprising, why does this method work in this case? I just tested this method on the linked problem, the same issue shown in the body of question comes up. What's the difference between these 2 problems?

@xzczd It works too in the linked problem, because point {.5,.5} is included in the definition of bmesh.

No, it doesn't. Please read the linked post carefully. There should be no peak at {0.5,0.5}.

@xzczd You must take a boundary point , for example {0,0} , for the additional DC-condition! Try U = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == f + NeumannValue[g, True] , DirichletCondition[u[x, y] == 1, x == 0. && y == 0. ] }, u, {x, y} \[Element] mesh] Plot3D[U[x, y] , Element[{x, y}, mesh]] which gives a smooth plot.

The issue is the same. This time the peak moves to {0,0}. Once again, please read the linked post carefully.

Can't believe it! Did you try my code with mesh = ToElementMesh[Rectangle[], "MaxCellMeasure" -> 0.001, "IncludePoints" -> {{0, 0}}, "MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];?

Yes, and the issue remains. (If it was that simple, we wouldn't have used such complicated methods to solve the problem. )

@xzczd, the only think I wanted to show with my comment is how to include a point inside the domain and how to address that point from a DirichletCondition. I am not claiming it solves any of these problems, I am just showing the technique.

@user21 I know. I'm just asking if you know why the method works in this case. You see, you're the one who knows FiniteElement best here :) .

@xzczd, hm I don't know. perhaps this has something to do with the fact that the disk is rationally symmetric. But this is really a shot in the blue. Sorry I don't know.