Approximating the laplacian and bilaplacian of a phase field

Dear NGSolve community,

recently, I got reason to numerically approximate high order derivatives (like the bilaplacian) of a phase field function.
To my surprise, I got into some trouble doing so and hope for being pointed out
a major fault in my method of using NGSolve or directed to a better method. The code I have been using to debug the issue is attached.

I noticed that the laplacian - and accordingly the bilaplacian - approximates show oscillatory behaviour. Sure, the phase field parameter I employ may be small and so I played around with making it larger or, on the other hand, refining the mesh while keeping it small.
Both approaches did not lead to a satisfactory result: The oscillatory behaviour still remains outside the
interfacial layer. To my understanding, computing gradients and divergences with NGSolve like in the attached code is based on differences of the phase field on a reference element. Am I mistaken here?

I am not very experienced in numerical issues so any recommendations will be appreciated.

All the best,
Philipp

Edit: I just substituted the phase field with the constant function 1 and even then the laplacian is not approximated correctly (Image attached), so there must be some error in using NGSolve I suppose.

if you interpolate into a p=1 finite element space, and then numerically differentiate 4 times, you cannot expect convergence. If you interpolate into a p=4 finite element spaces, you will see convergence.

well, the picture you are showing is zero up to 1e-11. The u.Set(1) does element-wise L2-projection, and you get some roundoff errors. Interestingly, by interpolating using dual shapes, the Laplacian will be exactly 0. You get that via

u.Set(1, dual=True)

Joachim

Thank you very much for the answer.

The hint with the polynomial degree is very helpful and has the desired effect in the setup I posted
here.
I recently tried to do computations with grid functions in 2d that are given as bitmaps.
The conversion into a grid function of order p=1 is no problem. Then I let NGSolve interpolate
into a higher order space by using the Set operation of a higher order grid function.
The approximations of Laplacian and Bilaplacian are then quite bad.
Is there an obvious flaw in this method?

I completely understand if you don't answer since this might not be NGSolve related. However,
any hints are of course appreciated.

All the best,
Philipp