In this paper, we present an original algorithm for the numerical simulation of tectonic movements and the related formation of geological faults. The approach is based on the use of Discrete Element Method, where the geological media is represented as an agglomeration of discrete particles which interact as elastic, visco-elastic, or elasto-plastic bodies. This approach naturally allows accounting for large deformations and discontinuities in the geological media; thus, allows simulation of faults formation. Implementation of the algorithm is based on the use of Graphical Processor Units. The Discrete Element Method requires a high number of floating point operations and logical operations per a single particle (degree of freedom) per time step; whereas the number of memory access operations is relatively low. Thus, the use of GPUs decreases the computational time substantially. Published under licence by IOP Publishing Ltd.