Since P1 finite elements are linear, the second derivative should be zero as you have computed it the first way. If you want a nonzero answer here, I think you will need to define your function on a higher-order FE space.

The second way you’ve programmed it, you are doing something different. First, you compute the gradient of u. Then you INTERPOLATE that gradient (which is piecewise constant in each component) onto new P1 FE functions. Finally, you compute the divergence of those functions to get your answer. I don’t think this is what you want to do…

If I understood your answer correctly, I have to convert my vector u from P1 to P2 finite element space, then compute the Laplacian operator using the first approach.

Code :

Vh1 (Th,P1);
Vh1 u ; //My vector
Vh2 (Th,P2);
Vh2 w;
w=v; //Conversion of the vector u from P1 to P2
Vh2 lap_U=dxx(w)+dyy(w);

If you actually need the value of the Laplacian, it would be better to solve the whole problem with P2 or higher-order elements and compute the Laplacian directly. Here, your code interpolates a bilinear function onto a biquadratic function, which could lead to spurious wiggles.