The advancement of scientific machine learning (ML) techniques has led to the development of methods for approximating solutions to nonlinear partial differential equations (PDE) with increased efficiency and accuracy. Automatic differentiation has played a pivotal role in this progress, enabling the creation of physics-informed neural networks (PINN) that integrate relevant physics into machine learning models. PINN have shown promise in approximating the solutions to the Navier-Stokes equations, overcoming limitations of traditional numerical discretization methods. However, challenges such as local minima and long training times persist, motivating the exploration of domain decomposition techniques to improve it. Previous domain decomposition models have introduced spatial and temporal domain decompositions but have yet to fully address issues of smoothness and regularity of global solution. In this study, we present a novel domain decomposition approach for PINN, termed domain discretized PINN (DD-PINN), which incorporates complementary loss functions, subdomain-specific transformer networks (TRF), and independent optimization within each subdomain. By enforcing continuity and differentiability through interface constraints and leveraging the Sobolev (H1) norm of the mean squared error (MSE), rather than the Euclidean norm (L2), DD-PINN enhances solution regularity and accuracy. The inclusion of TRF in each subdomain facilitates feature extraction and improves convergence rates, as demonstrated through simulations of two test problems: steady-state flow in a two-dimensional lid-driven cavity and the time-dependent cylinder wake. Numerical comparisons highlight the effectiveness of DD-PINN in preserving global solution regularity and accurately approximating complex phenomena, marking a significant advancement over previous domain decomposition methods within the PINN framework.