Hello Freefemers,
Testing the solver for Laplace equation with Neumann BC in 3d (in a cube). I observe a strange behavior: while the solution is perfect (error~10^-14) when using P2 FE, it is incorrect (error~10^0) when using P1 elements.
Of course, we expect a difference when using P1/P2 FE but not at at this order. So definitely there is an error with my code but I am unable to identify it. I have attached the code in the case anyone would like to have a look at it.
Thank you.
load “msh3”
int nn = 5;
mesh3 Th = cube(nn, nn, nn);
fespace Uh(Th, P2); // must be P2 for Neumann BC to work well
Uh u, ue, uh, f, err;
ue = -(x^2 + y^2 + z^2);
f = ue + 6; // f = ue - \Delta ue
macro grad(u) [dx(u),dy(u),dz(u)] //
solve neumann(u, uh) // solve u - Delta(u) = f, dn(u)=dn(ue)
= int3d(Th) ( u*uh + grad(u)'*grad(uh) )
- int3d(Th) ( f*uh )
- int2d(Th, 1,2,3,4,5,6) ( (N’*grad(ue))*uh ); // works well if P2
err = u - ue;
cout << "******* error = " << err[].linfty << endl;