We present a novel approach that integrates unfitted finite element methods and neural networks to approximate partial differential equations on complex geometries. Easy-to-generate background meshes (e.g., a simple Cartesian mesh) that cut the domain boundary (i.e., they do not conform to it) are used to build suitable trial and test finite element spaces. The method seeks a neural network that, when interpolated onto the trial space, minimises a discrete norm of the weak residual functional on the test space associated to the equation. As with unfitted finite elements, essential boundary conditions are weakly imposed by Nitsche’s method. The method is robust to variations in Nitsche coefficient values, and to small cut cells. We experimentally demonstrate the method’s effectiveness in solving both forward and inverse problems across various 2D and 3D complex geometries, including those defined by implicit level-set functions and explicit stereolithography meshes. For forward problems with smooth analytical solutions, the trained neural networks achieve several orders of magnitude smaller H 1 errors compared to their interpolation counterparts. These interpolations also maintain expected h - and p -convergence rates. Using the same amount of training points, the method is faster than standard PINNs (on both GPU and CPU architectures) while achieving similar or superior accuracy. Moreover, using a discrete dual norm of the residual (achieved by cut cell stabilisation) remarkably accelerates neural network training and further enhances robustness to the choice of Nitsche coefficient values. The experiments also show the method’s high accuracy and reliability in solving inverse problems, even with incomplete observations.