Summary: In fractured rocks, fluid flows mostly within a complex arrangement of fractures. Both the fracture network structure and its hydraulic properties are determined at first order by the broad range of fracture lengths and densities. To handle the observed wide variety of fracture properties and the lack of direct fracture visualization, we develop a general and efficient stochastic numerical model for discrete fracture networks (DFNs) in a three-dimensional (3D) computational domain. We present an original conforming mesh generation method addressing the penalizing configurations stemming from close fractures and acute angles between fracture intersections. Flows are subsequently computed by using a mixed hybrid finite element (MHFE) method. The lack of direct fracture knowledge is treated by Monte-Carlo simulations requiring simulations with a large number of networks with various characteristics. We analyze the complexity in size and in time for the computation of flow in 3D DFNs meshed with our method and compare with the complexities for 2D rectangular domains meshed with a regular grid. We find out that complexity in size is similar whereas complexity in time is slightly larger for DFNs than for 2D regular domains.