We present an algorithm for simulation of the Earth's crust tectonic movements and formation of the geological faults and near-fault damage zones. The algorithms are based on the Discrete Elements Method, and it is implemented using CUDA technology. We used to simulate faults formation due to different scenarios of tectonic movements. We considered the displacements with dipping angles varied from 30 to 90 degrees; i.e., up to vertical throw. For each scenario, we performed simulations for some statistical realizations. To characterize the simulated faults and damage zones, we consider the strains distribution and apply data clustering and Karhunen-Loeve analysis to distinguish between different forms of the fault zones. In particular, clustering analysis shows that displacements with high and low dip angles form completely different geological structures. Nearly vertical displacements, high dip angles, form wide V-shaped deformation zones, whereas the at displacements cause narrow fault-cores with rapidly decreasing strains apart from the fault core. Results of the presented simulations can be used to estimate mechanical and seismic properties of rocks in the vicinity of the faults and applied further to construct models for seismic modeling and interpretation, hydrodynamical simulations, history of matching simulation, etc.