  We present a new finite element method, called $\phi$-FEM, to solve
numerically elliptic partial differential equations with natural (Neumann or
Robin) boundary conditions using simple computational grids, not fitted to the
boundary of the physical domain. The boundary data are taken into account using
a level-set function, which is a popular tool to deal with complicated or
evolving domains. Our approach belongs to the family of fictitious domain
methods (or immersed boundary methods) and is close to recent methods of
cutFEM/XFEM type. Contrary to the latter, $\phi$-FEM does not need any
non-standard numerical integration on cut mesh elements or on the actual
boundary, while assuring the optimal convergence orders with finite elements of
any degree and providing reasonably well conditioned discrete problems. In the
first version of $\phi$-FEM, only essential (Dirichlet) boundary conditions was
considered. Here, to deal with natural boundary conditions, we introduce the
gradient of the primary solution as an auxiliary variable. This is done only on
the mesh cells cut by the boundary, so that the size of the numerical system is
only slightly increased. We prove theoretically the optimal convergence of our
scheme and a bound on the discrete problem conditioning, independent of the
mesh cuts. The numerical experiments confirm these results.
