Rock joints play a crucial role in the mechanical behavior of a rock mass and in the effective design of rock slopes and underground excavations. However, the shear behavior of rock joints is scale-dependent, and based on the previous studies, the in situ shear testing of rock joints both physically and experimentally is a huge constraint. In this paper, the shear behavior of natural jointed rock samples is captured through direct shear tests both in experimental and numerical modeling approaches. A discrete element modeling (DEM) code was developed in Particle Flow Code (PFC2D) where the intact rock is calibrated by a flat-jointed bonded particle assembly with a dense packing of circular particles bonded together. The results of the DEM approach are validated with the experimental direct shear test results. After the validation of the approach, a graphic processing unit (GPU)-based parallel computation algorithm was developed incorporating the direct shear test DEM code. The parallel computation of 22 GPU-accelerated nodes reduced the computation time significantly which helped to model rock joints of any size numerically. The scale effect on the shear behavior of rock joints is studied using this developed algorithm at different scales, and variation of the peak shear strength with respect to joint length is captured. Finally, it was concluded that the developed GPU-based parallel computation DEM code can be a useful approach for studying the scale effect on the shear behavior of rock joints without any size and computational time constraints.