I’ve seen a few posts in the past asking about implementing Neumann boundary conditions, especially homogeneous neumann boundary conditions in FreeFEM++.
I think it would be helpful to explicitly state somewhere that, since Natural boundary conditions vanish in your variational formulation, FreeFEM++ implicitly assumes any boundaries you don’t label are neumann.
Additionally, it might also help to describe how to implement Robin boundary conditions, where one implements it like:
If Th is a 3D mesh
\partial u / \partial n + g(x,y,z) = 0 on Gamma_1
You just put in