Abstarct: Connected discontinuities are known as main paths to flow fluids through the rock media. One of the most important modelling methods is the Discrete Fracture Network (DFN), in which the discontinuities and their connectivity patterns are considered as the only route for the flow, and the rock matrix is assumed impermeable. The geometry of DFN is established on partial differential equations, thereby statistical simulation of the fractures may resolve uncertainty problems to a large extent and represent more realistic models. These simulations can be even more precise using three-dimensional DFNs. The fractures are often modeled planar with a special – statistical position in the DFN simulation. The classic analytical relations have been developed to calculate the fluid flow in one-dimensional structures. However, these relations must necessarily be generalized to two-dimensional structures in order to simulate three-dimensional DFNs; therefore, the numerical discretization methods could be functioned. These methods reduce two-dimensional regions to linear elements and rebuild a new system of equations. One of them have been used from the recent half century due to its numerous advantages is well-known as the Finite Element Method (FEM). Ussing the FEM in discretizing the problem of the fluid flow in the fractured rocks result in a large sparse matrix called transmissevity matrix, the flow rate and the hydraulic head gradient vectors. Furthermore, particular methods are required to determine the model response. The solving methods of large sparse matrices are quite a lot various and classified to two main categories; direct and iterative methods. The direct methods are straightforward but they are not efficient enough to tackle all matrices. The most important iterative methods are Krylov Subspaces. Nevertheless, these methods have not well regarded in the rock engineering. In this research, the geometrical frxamework of the fractured rock medium is generated by the three DFN. In the sake for simulating the problem as real as possible, the effect of in-situ stresses on the fluid flow is considered, which leads to a one-way hydro-mechanic coupling and the calculations are implemented using the computer code, FlowSHUT3D. Validation of the results are performed by two different methods and sensitivity analyses are presented to comprehend the effect of key parameters on the fluid flow in the fractured rocks.