Numerical solvers for partial differential equations face a persistent trade-off between computational efficiency and accuracy, particularly for multiscale and time-dependent systems. Neural operators (NOs) have emerged as promising surrogate models due to their strong generalization capabilities and fast inference; however, their standalone deployment is limited by the need for large volumes of high-fidelity training data and reduced accuracy over large-scale domains. To address these challenges, we propose a hybrid domain-decomposition-based solver that couples physics-informed neural operators with the finite element (FE) method. In the proposed FE–NO framework, computationally expensive subdomains are delegated to a pretrained NO, while the remaining regions are solved using FE solver. Spatial coupling is achieved via the Schwarz alternating method, employing Dirichlet–Dirichlet conditions for overlapping interfaces and Neumann–Dirichlet conditions for non-overlapping interfaces. To ensure temporal coupling, we embed the Newmark–beta time-integration scheme directly into the NO, resulting in a time-marching NO that enables fully spatiotemporal FE–NO integration. Furthermore, by incorporating PointNet into the NO architecture, the resulting time-marching Point-NO can operate on arbitrarily shaped domains, allowing seamless replacement of irregular, complex, or highly refined subdomains within the hybrid framework. The proposed method is validated on a range of solid mechanics problems, including static linear elasticity, quasi-static hyperelasticity, and elastodynamics. In the hyperelastic setting, the FE–NO solver achieves a 20% reduction in computational cost compared to a conventional FE–FE solver, while maintaining errors below 5%. For dynamic simulations, the non-overlapping FE–NO formulation converges within only three inner iterations per time step, compared to nine iterations required by the overlapping strategy. Importantly, in the non-overlapping case, errors remain bounded across all 50 time steps, demonstrating negligible accumulation despite the auto-regressive time-marching scheme. Finally, we demonstrate the applicability of the framework to a realistic, multiscale, and highly nonlinear problem: the mechanical response of lung tissue under blast and blunt impact loading. This hybrid solver reduces wall-clock simulation time from days to hours, highlighting its potential for large-scale, high-fidelity computational mechanics application.