Direct Numerical Simulation (DNS) is the most accurate method in computational fluid dynamics (CFD) as it directly resolves all turbulence scales, eliminating the need for turbulent modeling approximations. Despite its accuracy, DNS is impractical for many applications due to the extremely fine spatial and temporal resolution required. However, with the recent development of high-order numerical methods such as Continuous Galerkin (CG) and Discontinuous Galerkin (DG) spectral element methods, the grid resolution required is significantly reduced compared to lower-order methods, improving the feasibility of DNS in certain scenarios. These methods achieve high-order accuracy by utilizing high-degree polynomial basis functions. Nevertheless, they face significant limitations in handling numerical solutions with sharp gradients and discontinuities, such as those seen in turbulent flows and shock waves. This challenge is due to the inherent difficulties in capturing these features with polynomial functions. To address these challenges, I will introduce nodal integral methods, which combine polynomial and non-polynomial basis functions derived from semi-analytical solutions, providing better coupling of physical processes and potentially enabling more feasible DNS simulations.