How to solve Laplace equation with Neumann boundary condition?

Hello everyone,
I am using to Freefem to solve a very simple equation: Poisson equation with Neumann boundary condition. I have read the document, but it just said about the Dirichlet example! I don’t know how to put the Neumann boundary condition into the code! The document says that the neumann boundary condition will appear in the bilinear form but what if the Neumann boundary is 0, then the bilinear form is similar to the one with Dirichlet condition. I am very confused!
Thanks for any comments!

The Neumann bcs enter the variational form as an additional integral operator along the boundary at which they are imposed.

… + int1d(Th, ID) (v_test*du/dn),

where du/dn is the known Neumann boundary data you want to impose.

Thanks for your comments! But what if du/dn equals 0? Do we still need to put the int1d there? In addition, what about the uniqueness of the poisson equation? From what I learned, the solution u is unique up to a constant! So, how we fix the constant in the code?
Thank you very much!

First the laplace bvp with only Neumann data is not uniquely solvable (ill-posed). In particular you can verify that any constant function satisfies your problem, so in essence you do not need to solve anything; just pick one constant (although I do not see why you would like to do that). In general if you have homogeneous Neumann data the boundary integral term just vanishes, and hence the Neumann BC is again correctly imposed.

Thanks a lot! Please help me to answer one more question! You said that you can pick any constant. For example, you solve this equation -Delta u = f and du/dn = g, if u is a solution then u + c is also a solution! if f is 0, then you can pick any constant but what if f is not 0? I think that I have seen the code of solving the Stokes equation in Freefem, it has a similar issue with this neumann problem. There, they put 1e-06pq to fix the constant term, so I just wonder if we do the same thing for the neumann problem?
Thank you for your answers!

As you say the Poisson problem is also not uniquelly solvable, when employing only Neumann data. You also need the compatibility condition

int2d(Th)(f)+int1d(Th,Boundary_ID)(g_N) = 0

The constant can be fixed by adding constraints; for instance you can incorporate the constraint that the average of your solution is vanishing, using for instance a lagrange multiplier approach to impose the constraint int2d(Th)(u)=0 or other regularization terms.

Regards

Thank you very much for reminding me about the compatibility condition! In analysis, we can impose the average of u = 0 but my question is how we put this impose into the codes? Should we add epsilon*(u,v) into the bilinear form? I have seen this technique from Stokes equation! The pressure term is unique up to a constant! Hence, they add epsilon*(p,q) to stable the scheme! However, I am not sure that we can do the same here!
Thank you very much for your helping me! I am really appreciate that!

Please have a look at this PETSc example or this sequential example.

Thank you very much! The example is useful!