Eikonal tomography, or travel time inversion, has been one of the primary seismological tools for decades and has been used to understand Earth’s properties and dynamic processes. At the heart of the inversion process is the need for an accurate, and preferably flexible, eikonal solver to compute the travel time field. Most of the conventional eikonal solvers, however, suffer from first-order convergence errors and difficulties in dealing with irregular computational grids. Physics-informed neural networks (PINNs) have been introduced to tackle these problems and have shown great success in addressing those challenges. Nevertheless, these approaches still suffer from slow convergence and unstable training dynamics due to the multi-term nature of the loss function. To improve on this, we propose a new formulation for the isotropic eikonal equation, which imposes boundary conditions as hard constraints. We employ the theory of functional connections to the eikonal tomography problem, which allows for the utilization of a single loss term for training the PINN model. Through rigorous numerical tests, its efficiency, stability, and flexibility in tackling a variety of cases, including topography-dependent and 3D models, are attested, thus providing an efficient and stable PINN-based eikonal tomography.